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Abstract 


This thesis models analytically and numerically a moving boundary prob- 
lem originating in a small particle accelerator known as a neutron tube when 
a positive-ion plasma is exposed to a large accelerating electric field. 

The coupled non-linear system of equations describing the flow of charge 
within the acceleration region, and the consequent plasma boundary loca- 
tion, are modelled in both one and two dimensions. The one-dimensional 
study provides an analytic solution to the steady-state problem in planar 
geometry that arises with given, steady boundary conditions, and goes on to 
develop a successful numerical method for the determination of a numerical 
solution to the same problem. The one-dimensional study continues with 
an analytic analysis of the time-dependent system, concluding with the de- 
velopment of numerical methods devised to analyse the full time-dependent 
problem. The two-dimensional work builds on the successful one-dimensional 
time-dependent numerical method, using space-charge limitation to aid in the 
determination of boundary location techniques, thus developing a model for 
the steady-state two-dimensional case. The method initially utilises a fixed 
boundary, iterative method as a vehicle for calculating the initial location of 
the boundary, and continues, having adjusted the boundary (and hence the 


solution region), until a termination criterion is satisfied. 
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Chapter 1 


Introduction 


1.1 Problem Background 


Neutron generators have been used as exploratory tools for a number of years, 
with the highly sensitive specific technique of Neutron Activation Analysis 
(NAA) being used to determine the elemental content of materials to a high 
degree of accuracy (See [39]). The method of NAA involves the illumination 
of the required sample by a pulse (or pulses) of neutrons resulting in the 
release of a characteristic signature from the sample, in the form of other 
types of radiation. 

During the neutron illumination process, the sample becomes activated 
and subsequently emits secondary gamma radiation in two phases, as it de- 
activates. The first phase is a rapid de-excitation of the activated nuclei via 
gamma decay where the characteristic gamma frequencies are particular to 
the element originally excited. The second phase is a radioactive decay pro- 
cess from unstable isotopes formed from the neutron illumination process, 


the decay rates being strongly dependent upon the isotopes formed (half 
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lives ranging from 10~° to 10” years). It is usually this second phase that is 
analysed to determine information about the sample constitution. 

For many applications of NAA, the neutron source is a large, laboratory- 
scale piece of apparatus. However, in some applications (such as oil well 
logging), such a large scale piece of equipment cannot be used, and so the 
neutron source must be taken to the location of the sample. In this case, 
the neutron source is a much smaller device and is a neutron generator, 
comprising a power supply and sealed neutron tube. 

The problem analysed in this thesis originates within a neutron tube 
that has an entirely different application, where it is used in nuclear weapon 
initiation. Briefly, the short, very intense pulse of neutrons created by the 
neutron tube is used to initiate a fission chain reaction within a nuclear 


weapon! [40]. 


1.1.1 Neutron Tubes 


1.1.1.1 Fusion Reactions 


2 


In order to generate neutrons within a neutron tube, fusion“ reactions are 


employed. By utilising fusion, the creation of neutrons can be precisely timed 


‘A more detailed description of this application cannot be given due to classification 


issues. 


?If fission reactions were to be used as a neutron source, a donor reactant would be 
chosen whose spontaneous fission products would be introduced to another fissile material. 
Subsequently, a further fission reaction would take place and one of the products of this 
second fission reaction would include neutrons with a specific energy. Creating neutrons 
by fission is less controllable by virtue of the spontaneous decay of the fission reactants, 


and on the chemical mixing of such materials. 
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and controlled by initiating the reaction only when required. 

To create the fusion reaction within the tube, an ion beam consisting of 
one of the fusion reactants is accelerated to a relatively high energy towards 
a stationary target that is infused with the other reactant. The energy to 
which the incident particles are accelerated is chosen depending upon the 
fusion reaction cross-section, and is dependent upon the particular reaction 
concerned. The reaction cross-section relates to nucleon - nucleon scattering, 
and is a measure of the probability of a particular reaction taking place (see 
[21], and Appendix A). It effectively gives an indication of the apparent size 
of the target nucleon for a given reaction and given energy. The reaction 
cross-section has the dimensions of area and is usually measured in barns 
where 1 barn is 100 fm? or 10~73m?. Sometimes, the reaction cross-section 
is given in a centre of mass frame of reference. In this reference frame, the 
reaction probability is independent of the total momentum of the reacting 
particles relative to the laboratory. The D-T cross-section referred to in this 
thesis is measured in the laboratory frame of reference, where the target 
particle is stationary relative to the laboratory. It is scaled from the generic 
centre-of-mass frame of reference cross-section (for the two given reactants) 
according to the reacting particle masses (Appendix A). 

Most fusion reactions capable of creating neutrons have negligible reaction 
cross-sections (i.e. they do not generate an appreciable number of neutrons 
for a given energy) unless the incident ion is accelerated to such high energies 
that the process is virtually impossible to achieve in a small neutron tube. 
This is since the large electrostatic potential used to accelerate ions must 


be held-off across a small distance of the order of mm. The term hold-off 
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refers to the high voltage electrodes within the neutron tube being separated 
by an insulating wall. The insulator used is usually a ceramic that not only 
separates the electrodes and hence resists electrostatic breakdown (where the 
potential energy of electrons at the cathode is so great that they are able to 
jump the gap between the electrodes either directly across the vacuum gap, 
or by skipping along the ceramic wall), but it also maintains a seal that 
resists permeation from gases within the atmosphere, thereby maintaining 
the vacuum within the sealed tube. 

To utilise low cross-section reactions, not only would the device have to be 
relatively large (scaled by, say, a factor of ten) and cumbersome, but the tube 
power supply would need to generate the huge potential differences necessary 
to drive the reaction; this is very difficult to achieve. Therefore, it is usual 
for a neutron tube to utilise the deuterium-tritium fusion reaction (or D-T 
reaction), whereby deuterium ions are accelerated on to a stationary (relative 
to the laboratory) tritiated target. This reaction has a relatively high peak 
cross-section, which occurs when the incident deuteron has a kinetic energy 
of ~110KeV. It is noted that the converse T-D (tritium-deuterium) reaction, 
whereby tritium ions are accelerated onto a stationary deuterium target, 
has a peak cross-section when the incident triton has a kinetic energy of 
~165KV (see Appendix A). The laboratory frame reaction cross-sections for 
these reactions are shown in Figure A.1. 

For the D-T fusion reaction to take place with the maximum probabil- 
ity, deuterium ions must be accelerated to an energy of ~ 110KV before 
striking the tritium impregnated target. Upon striking the target, some of 


the incident deuteron particles fuse with the tritium within the target (the 


ral 


proportion of fused incident ions being dependent upon the reaction cross- 
section) consequently releasing neutrons with 14.1MeV in addition to 3.5MeV 
alpha particles. Those deuterium ions that do not fuse with tritium lose en- 
ergy as they penetrate the target, thereby losing some of their ability to fuse 


(by effectively sliding down the cross-section curve) at the next interaction. 


1.1.1.2 Tube Operation and Plasma Boundary Formation 


Geometrically, the tube consists of a sealed envelope, across which the ~ 
110KV electrostatic potential difference is held off by a ceramic wall. A source 
of deuterium ions is situated at the anode end of the tube whilst the tritiated 
target is at the cathode end (see Figure 1.1). Upon operation, a plasma of 
deuterium ions expands into the acceleration gap acting as a conductive 
‘gas’ which is generally impenetrable to the tube main accelerating field (due 
to its conductivity). At about the same time, the ~ 110KV accelerating 
voltage is applied across the tube and ions begin streaming away from the 
plasma-vacuum interface (the plasma boundary) formed by the expanding 
plasma. Shaped electrodes within the tube act as ion lenses focussing the 
ion beam onto the tritiated target where the fusion reaction takes place. 

There are two distinct timescales that exist within the operation of the tube. 
The timescale over which the source operates and the accelerating voltage is 
applied, is generally several orders of magnitude greater than the timescale 
over which the plasma boundary settles to a stable position for a given set of 
conditions. The accelerating voltage is ramped up in a period of jus, whereas 
for a given impulse, the boundary settling process occurs in a period of ns. 


Since the tube operates in a dynamic manner, whereby the source and 
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Figure 1.1: Neutron tube schematic. 


accelerating voltage are switched on, reach relatively stable conditions and 
are then switched off, then due to the relatively long time scales of the source 
drive and acceleration voltage pulses (in comparison to the timescale of the 
boundary settling process), it is believed that the boundary moves in accor- 
dance with the applied acceleration voltage, eventually settling to a stable 
location during the relatively stable phase of tube operation. The boundary 
movement process takes place in the following way. 

Ions reaching the plasma-vacuum boundary at a specific rate cause the 


boundary to bulge into the vacuum (see Figures 1.2 and 1.3 where the ion 
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Figure 1.2: Schematic showing ions reaching the plasma - vacuum boundary prior 
to reaching the equilibrium state. In this state, the steady (relative to the timescale 
of the settling process) current density, Jin(x,y), of ions reaching the boundary 
from below is higher than the current density, Jout(x,y,t), of ions leaving the 
boundary from above. The horizontal lines are simulated level curves of electric 


scalar potential. Note that the region is symmetric about the vertical axis. 


arrival rate and departure rates are represented by the current densities 
Jin(x,y) and Jout(x,y,t) respectively), thereby concentrating the electric 
field within the region (since the plasma acts as a Dirichlet boundary condi- 
tion for the electric scalar potential within the tube main gap). The increased 
electric field (as a consequence of the reduction in electrode separation) causes 
ions to be accelerated away from the boundary more rapidly than they ar- 
rive there, and consequently the boundary recedes until the electric field at 
the boundary becomes zero [8] (see Figure 1.4). In this equilibrium state, 


ions leave the boundary at the same rate as they arrive there, and the the 
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Figure 1.3: Schematic showing ions reaching the plasma - vacuum boundary prior 
to reaching the equilibrium state. In this state, the current density, Jin(x,y), 
of ions reaching the boundary from below is now lower than the current density, 
Jout(x, y,t), of the ions leaving the boundary from above, since the electric field 
strength at the boundary has been increased by it bulging into the region, hence 
causing ions to be extracted more rapidly. The plasma - vacuum boundary now 


begins to recede. Simulated level curves of electric scalar potential are also shown. 


plasma-vacuum boundary position has stabilised. 

The determination of the position of the plasma-vacuum boundary during 
the settling process is a moving boundary problem, in which the boundary 
converges to an equilibrium position over a short period of time. On the 
timescale of the overall tube operation, the determination of the plasma- 


vacuum boundary position for a given source current density? and accelera- 


3The source current density profile is determined by the source metallurgy and geometry 


in addition to currents and voltages applied to it. 
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Figure 1.4: Schematic showing ions reaching the plasma - vacuum boundary at 
the equilibrium state. In this state, the current density, Jin(x,y) of ions reaching 
the boundary from below is equal to the (now stable) current density, Jout(x,y), 
of ions leaving the boundary from above. Simulated level curves of electric scalar 


potential are also shown. 


tion potential difference is also a free boundary problem. It is the purpose 
of this analysis to determine the free boundary positions for problems with 
given steady state source current densities and acceleration voltages using 
mathematical techniques tailored to the problem. In addition, the full time 
dependent moving boundary problem occurring at the shorter timescale is 


also explored. 


1.1.2. Previous Work 


The problem of determining the position of the plasma boundary within 


an ion accelerator setup is one that has attracted a good deal of attention, 


26 


particularly in the 1960s and 1970s. As a consequence, I have had to deal 
selectively with the references in an attempt to develop a coherent account of 
past work, without reviewing subject matter that does not directly influence 
the current work (although some of this experimental work is mentioned in 
passing). Here I present a general history of related work, but detail on some 
of the more relevant works is given in Chapter 2. 

Tonks and Langmuir [2] first theoretically analysed the behaviour of a 
plasma expanding into a region with an applied electric field by investigating 
low pressure plasma discharge within a discharge tube (a plasma is defined 
as a neutral gas consisting of electrons and positive ions). In this work, 
the ion number density distribution within such a plasma is determined by 
approximating solutions to specific cases of what they called the ‘complete 
plasma sheath equation’ (an integro-differential equation) in one dimension 
for different geometries (more general solutions of this equation were also 
considered by Self |7]). 

The analysis of Tonks, Langmuir and Self detailed specifically the be- 
haviour of an expanding plasma up to and including the plasma sheath? 
(or the plasma boundary), but their analysis was unable to adequately de- 
scribe the self-consistent electric field outside the plasma region (beyond the 
sheath). This was primarily due to their representation of the ion number 
density (the free charge density of ions divided by their charge) within the 
plasma being no longer valid outside the plasma region. More recently, in an 
independent work ([{11]), a similar approach (to Tonks and Langmuir) was 
adopted, but the ion number density within the plasma region was calculated 


4The region where electrons within the plasma are repelled by the external electric 


field. 
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from the ion equations of motion instead of using the approach of Tonks and 
Langmuir. This enabled the calculation of the correct potential distribution 
outside the plasma region, whilst also determining the position of the plasma 
sheath, although the ion number density within the plasma was not correct. 
This is true because the representation of ion acceleration outside the plasma 
region is correct, but inside the plasma region it is not correct. The discon- 
nected nature of the plasma region and ion acceleration region enables this to 
be the case. Electric scalar potential distributions for a specific application 
were then calculated in one dimension. As a way of determining the plasma 
sheath thickness and position a Maxwell-Boltzmann distribution ([18]) was 
used to represent the electron number density within the plasma (as in pre- 
vious works) and solutions were approximated by asymptotic methods. 
Other work detailing the behaviour of an expanding plasma originate from 
the experimental discovery by Tonks and Langmuir in 1929 ([|1]) of plasma 
oscillations (the ability of a neutral plasma to support wave motion and 
potentially shock like behaviour). Previous models (such as those mentioned 
above) of the plasma behaviour were unable to support such oscillations. The 
inclusion of the Vlasov equations ([41], [32]) (a hyperbolic system describing 
the behaviour of the ion number density in a plasma with a long range 
charge interaction) as a description of the ion charge density would enable 
such plasma oscillations to occur within the plasma model. Works such 
as Whealton et al [12] and Whealton, Bell et al [22] solve the non-linear 
Vlasov-Poisson system (a similar system to that originally proposed by Tonks 
and Langmuir but where the ion number density is determined from the 


Vlasov equation) and consequently trace ions through the plasma and into 
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the electrostatic field in a Lagrangian manner. 

Treating ions as particles and tracing them through calculated field so- 
lutions is a method commonly used for the determination of the ion charge 
density outside the plasma region in ion beam numerical methods (refer, for 
example, to [28] and [31] in addition to references mentioned above). There 
are several problems associated with this method of predicting the ion beam 
location, the most notable being the number of ion beams that must be used 
to accurately represent the charge density distribution within the solution 
region. Not only is accurate ion beam tracing very computationally inten- 
sive, but accurate methods of assigning the calculated space charge density 
to the field solver mesh must be used as part of the computation procedure 
([27)). 

Most field solvers used in ion beam tracing codes have used static mesh 
finite difference approximations to the field equations. However in 2004, 
Humphries [37] proposed a free surface finite element approximation to the 
plasma boundary problem, whilst maintaining the Lagrangian particle trac- 
ing procedure. In this work, the free surface position is determined by main- 
taining a constant ion current density at the surface (thereby placing the 
required additional boundary condition there). Nodes along the emission 
surface are iterated to a stable position by moving them at each iteration 
with the distance moved being determined from how far away from the re- 
quired constant current density the calculated emission current density is. 
For a more in depth history of work in this area (until 1983) refer to Lejeune 
[16] where a very comprehensive set of references is also listed. 


Our interest essentially follows on from [11], but we are not interested 
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in the plasma or plasma sheath per-se (as many previous works are). Con- 
sequently we concentrate on the boundary value problem that the plasma 
expanding into an electric field presents. Since we are not interested in the 
plasma, but merely the plasma boundary position, we do not include the 
electron number density in calculations and seek the plasma boundary po- 
sition as a consequence of other conditions. The known (from experiment) 
current density distribution at the plasma boundary is used (this is not nec- 
essarily constant in position) in conjunction with a calculated ion charge 
density (from the ion equation of motion), and extra Neumann boundary 
conditions (see [8]) at the plasma boundary. We calculate the electric scalar 
potential distribution in addition to the charge density and charged particle 
velocity distributions, treating the particle velocity as a conservative field. 
Also, the position of the plasma boundary position or the ion emission region 
shape (treating this boundary as the edge of the solution region) is calcu- 
lated both analytically and numerically in one dimension and numerically in 


two-dimensional planar geometry. 


We begin by discussing in detail, the formulations of Tonks and Langmuir 
and of the internal work (these are most relevant to the model adopted in this 
thesis), [11], in addition to some of the other work mentioned above. Subse- 
quently, we formulate a mathematical description of the problem investigated 
here, discussing methods employed to solve the problem. Our analysis be- 
gins with a solution of the simpler one-dimensional steady state problem, and 
follows with a full time dependent solution of the one-dimensional problem 


with a moving boundary. Following on from this, using information acquired 
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in the one-dimensional case, an approach to solving the more difficult two- 
dimensional case is given whereby the time dependency in the field equations 
is used in a pseudo manner as a device for obtaining convergence in the free 


boundary case. 


3l 


Chapter 2 


Mathematical Formulation 


The situation outlined in Chapter 1 corresponds to a system of coupled prob- 
lems describing the various phenomena occurring within a neutron tube dur- 
ing normal operation. These phenomena, as previously indicated, can be 


broken down into four main areas and are listed as follows. 


1. Source Operation and Plasma Expansion 


The tube ion source produces a plasma consisting of both positively 
charged deuterium ions and electrons. After release from the ion source, 
the plasma expands into a field free region (termed the plasma cup) 
until it meets the high voltage accelerating field by passing through an 


aperture. 


2. Plasma - Vacuum boundary 


As the neutral plasma reaches the high voltage accelerating field re- 
gion of the neutron tube, a sheath forms at the boundary between the 
neutral plasma and acceleration region. Within the sheath, the den- 


sity of electrons gradually reduces with distance into the acceleration 
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region until all electrons are repelled by the electric field. We define 
this location as the plasma - vacuum boundary, and it corresponds to 
a boundary shape whereby the flux of ions reaching it from the source 
is matched by the flux of ions leaving it towards the target. It is also 
defined by a zero field condition (see |8]) since if a non-zero electric 
field were to exist on the acceleration field side of the boundary, then 
ions will be extracted from the boundary more rapidly than they arrive 
there from the plasma side and the boundary would adjust until the 
balance is restored. The shape of the boundary is dependent upon the 
deuterium ion current density across it, the potential difference across 
it, and the geometry of the acceleration region, in addition to the zero 


field condition at the boundary. 


3. Acceleration Region 


Once ions reach the plasma boundary, they are extracted and acceler- 
ated towards the target. The electric field within the acceleration region 
is dependent upon the charge density within it, the potential difference 
applied across it, and the shape of the plasma boundary formed at the 
anode. The ion beam shape is also dependent upon the electric field 


and hence its own charge density. 


4. Target 


After being accelerated across the neutron tube, the (now) high energy 
ions strike the tritiated target undergoing fusion, and release 14.1 MeV 


neutrons isotropically. 


TA plasma by definition is a field free region. 
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For reasons that will become more apparent in the following sections, it is 
possible to decouple the behaviour of ions within the plasma from ions within 
the acceleration region. Since we are not specifically concerned with the 
modelling of the plasma itself, it is only the acceleration region and position 
of the plasma boundary that is of concern in this thesis. However, two plasma 
models that lead logically on to the boundary value problem considered here, 


are now discussed. 


2.1 Tonks and Langmuir 


As mentioned previously in Section 1.1.2, Tonks and Langmuir ([2]) consid- 
ered the ion density within an expanding plasma and formulated what they 
termed the complete plasma sheath equation, for a one-dimensional plasma 
expanding from the origin into a region with a potential difference applied 
across it (see Figure 2.1). This equation was constructed from the Maxwell 
continuity equation (or Gauss’ Law) relating the divergence of the displace- 
ment current from a region to the charge density enclosed within that region. 


In one dimension this can be written as 
_ | 
es ee NSN ly (2.1) 
0 


where ¢(x) is the voltage (electric scalar potential), g the electron charge, € 
a scaling parameter called the permittivity of free space, and N; and N, the 
ion and electron number densities within the plasma, respectively (here it is 
assumed that ions are singly charged). 

Equation (2.1) is written in this form (the charge density as the product 


of particle charge and number density) to allow the use of the Maxwell- 
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Figure 2.1: The region upon which the complete plasma sheath equation is solved. 
The Dirichlet voltage conditions $(0) = ¢9 and ¢(r) = ¢1 are applied at the region 


boundaries. Since ions are positively charged, ¢o > ¢1. 
Boltzmann distribution ([{18]), 
f(T, €) = Noe, (2.2) 


to describe the electron number density within the plasma. The Maxwell- 
Boltzmann distribution? function f(T, ¢) gives the mean number of particles 


in a state of energy ¢€ for a system of particles at absolute temperature T (k 

*This distribution is applicable when the particles are completely distinguishable from 
one another. This in the absence of quantum effects which become prevalent when the 
particle density is such that the mean distance between the particles is of the order of 
their De Broglie wavelength. When quantum effects become important, their number 
distribution for a given energy depends upon whether the particles are bosons (Bose- 
Einstein distribution with integer spin quanta) whose wave-functions can superpose, or 
fermions (Fermi-Dirac distribution with half integer spin quanta) whose wave-functions 


cannot superpose. 


39 


being the Boltzmann constant). This is the most likely number of particles 
with that energy and that temperature, and is determined by analysing the 
number of ways of arranging particles in a particular state within the system, 
and choosing the state with the greatest number of arrangements possible 


(and hence the most likely state). Setting 


= Noe®r (2:3) 


from (2.2), where e(¢) = —q@, and where Np is the total number of electrons 
within the plasma, gives an expression for the electron number density as a 
function of voltage in the Poisson equation (2.1). 

If we assume that ions are created (atoms are ionised) throughout the 
expanding plasma (an assumption that would be valid for a discharge ion 
source), an expression for the total ion number density at a given position x 
within the region can be determined from the ion number density per unit 
length n(x) at that position. This is related to an ion generation rate per 


unit volume G(x) at x by the relationship (in one dimension) 


where v(x) is the ion speed at x. The total ion number density N;(x) at x 


gia * G(z) 7 
No) = fo Se ae, (2.4 


with the ion speed at z being determined by conserving energy from 


can then be written 


smv(z)” = 4 (9(z) — o(2)), (2.5) 
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where $(z) is the voltage (potential energy per unit charge at z) and m is 


the ion mass; rearranging (2.5) gives the ion speed 


v(z) = 4/> (o(z) - 6(2))?. (2.6) 


Substituting (2.3), (2.4), and (2.6) into (2.1) gives the complete plasma 


sheath equation 
Po _ 4 &) *_ G2) ange? a7 
da? Al DD tecane ny 


This equation describes the voltage profile, ¢(a), within a plasma in addi- 


tion to the ion and electron number densities as a consequence of the applied 
electric field. The interface between the plasma and vacuum region (Figure 
2.1) is determined as the location where the electron number density density, 
(2.3), is less than some given parameter 6, where 6 < No. At this location, 
effectively all electrons have been repelled by the external electric field (this 
point is not distinct as the electron number density decays to zero exponen- 
tially as the voltage decreases with the dependent variable x). Beyond the 
plasma-vacuum interface, no further ions are generated and the integral in 


(2.7) becomes 


where C' is a constant depending on the position s of the interface. Addi- 


tionally, since the electron number density in the region beyond the interface 
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is effectively zero (that is, less than 6), the incorrect implication is that the 
ion charge density outside the plasma is constant. Furthermore, Tonks and 
Langmuir sought to model a discharge ion source where ions are generated 
throughout the plasma region. The ionisation region is localised in modern 
neutron tube designs and so such a model is inappropriate. We therefore 
conclude that (2.7) is not adequate for the problem at hand. To correctly 
describe the ion beam behaviour after ions have been stripped away from the 
plasma surface, a different model of the ion density in this region must be 
employed, and this is most readily achieved using the equations of motion 


for the ions as they are accelerated by the external electric field. 


2.2 Internal Model 


The internal model, [11], investigates the one-dimensional expansion of a 
plasma within the region x € (—oo, 0] with the applied Dirichlet conditions 
@(—oo) = 0 and (0) = do, where ¢(x) is the voltage at x, as before. In 
order to describe the voltage distribution within this region, the same form 
of the Maxwell equation (2.1) used by Tonks and Langmuir is employed. 
However, instead of modelling the generation of ions throughout the plasma, 
it is assumed that ions are generated only at the left hand end of the region 
(at x = —oo) and their number density is calculated from their speed at 
x, v(x). It is noted here, that whilst ¢(0) = ¢9 and v(—co) = vp seems 
contradictory, it is the notation used in the original paper. 


Equating the ion kinetic and potential energy (from the applied electric 
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field) at x we have 


5 (v(2)? — vj) = —4¢(6(x) — ¢(—00)) 
= —q¢(z), 
v(x) = 4/ ve — “1, (2.8) 


where m, q and v(x) are as previously defined, and where vg = v(—oo). 
Denoting N;(a) as the ion number density, as above, then the ion flux I(x) 


at x is given by? 


In a steady state of ion flow 
dV 
ee 
dx ; 
which implies that [is independent of x, or that 
T(z) = Novo 


= const, (2.9) 


where No = N;(—oo). This gives the expression 


Novo 
v(x) 


N,(2) = 2, (2.10) 


from (2.8) and (2.9) for the ion number density. Since ions are created 
by ionising neutral atoms, then by conserving charge we can also say that 
N-(—oo) = No. By assuming the Maxwell-Boltzmann distribution (2.2) for 


3The one-dimensional definition of flux density. 
4From the continuity equation V -T = ee where the time derivative is zero. 
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the electron number density (with N.(T,¢) = f(T, ¢(@)) as previously de- 


fined), the Poisson equation, (2.1), can be written 


b> Noa f (2 2q_\7? 
$e = Ait |e) («3 *46) \ (2.11) 


E0 
from (2.2) and (2.10). Furthermore, by defining the dimensionless parameters 
mug 
oi 


ae 
RE’ 


where 


is the Debye shielding length ([41], [32]), (2.11) can be re-written in the scale 


eae (1-2) (2.12) 


dz? a 


invariant form 


NI 


Solutions to (2.12) were calculated for a typical neutron tube application. 
These solutions gave a good indication of the location of the plasma boundary 
based upon the electron number density and the calculated self-consistent 
accelerating potential. Results from this model indicate that ions within the 
plasma are shielded from the external accelerating field by the presence of 
neutralising electrons. The shielding shows that within the plasma region, 
ion trajectories are completely unaffected by the external field, and only 
when the electron number density begins to drop within the thin plasma 
sheath region do the ions begin to feel the external field. This leads to the 


concept (as it was termed) of the plasma “freezing in” the ion trajectories 
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due to its self-shielding effect, or that the plasma and acceleration regions 
are effectively uncoupled except for the very thin plasma sheath region. 
This model leads nicely on to the present work, where we assume that 


the plasma and acceleration regions are uncoupled. 


2.3. Current Model 


A representative neutron tube acceleration region is shown in Figure 2.2, 
with the ion source located at the bottom and the target at the top (both in 
black). The plasma boundary (marked) is situated immediately above the 
aperture in the ion source plasma cup. A representation of the ion beam is 
shown by the two curved lines separating the zero charge density and non- 
zero charge density regions, stretching from the ion source cup to the target. 
The insulators provide resistance to gas permeation, in addition to holding 
off the voltage between the source and target electrodes; they have a different 
dielectric permittivity to the vacuum they enclose. 

This schematic is representative of either a two-dimensional planar re- 
gion, or as a cross-section through a three-dimensional axially symmetric 
arrangement. In both cases, a vertical line of symmetry (marked) separates 
identical parts of the region on either side of it. In the axially symmetric 
case, the vertical line represents an axis of rotational symmetry thereby re- 
ducing the dimension of the problem (since only axial and radial coordinates 
need be represented); in the planar case, the vertical line mirrors the solution 
on either side of it. With this in mind, if it is assumed that the tube insu- 


lators are a sufficient distance away®, then depending upon the numerical 


By sufficient distance, we mean far enough away that they do not have any influence 


Al 


algorithm chosen, it is only the evacuated region (zero and non-zero charge 
density) immediately above the plasma boundary that need be modelled in 
order to fully determine the ion beam voltage and velocity distributions. The 
ion beam profile and consequent boundaries forming the edge of the ion beam 


should naturally emerge as part of the solution. By observing the vertical 
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Figure 2.2: Tube schematic showing the ion beam (marked ’non-zero charge den- 
sity region’), the charge free region (marked ’zero charge density region’), and the 


plasma boundary situated directly above the source aperture. 


line of symmetry at the centre of the ion beam, the region of interest shown 
in Figure 2.2 can be represented as Figure 2.3 (this will be referred to in the 
ensuing formulation), where the region boundary is split into seven distinct 


sections; these are labelled {S;}. 


on the electric fields in the vicinity of the ion beam. 
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Figure 2.3: Region of interest, R. The curved line is a representation of the edge 


of the ion beam, and separates the zero and non-zero charge density regions. 


2.3.1 Field Equations 


The following set of equations describe the ion beam and voltage profile in 
the evacuated region shown in Figure 2.3. As described, due to the nature 
of the problem, the system is dynamic with the flow of charged particles 
being governed by the current density at the plasma boundary and by the 
applied potential difference (the shape of the plasma boundary itself also 
being dynamic); consequently, the electric field within the region is also dy- 
namic. Therefore, the set of equations describing the problem is initially 
time-dependent, later becoming time-independent upon settling to the equi- 


librium state. 


43 


2.3.1.1 Time-dependent, Unsteady Equations 


The full set of equations describing the problem in the time-dependent regime 


are 

Ov 

mi (Vy Oh = q(B+¥ xB), (2.13a) 
_ _ Op 
V -(pv) = OE (2.13b) 
VES, (2.13c) 
Eo 
V-B=0, (2.13d) 
OA 

and, B=VxA. (2.13f) 


Here the constants m and q are the ion mass and charge respectively (as in 
Tonks and Langmuir [2], and in the internal model |11]), the vector v(x, t) 
is the velocity vector field generated by the passage of particles across the 
region under consideration (x represents the vector of linearly independent 
coordinates appropriate for a given coordinate system), and the scalar field 
~(x,t) is the voltage within the region (as above). The quantity p(x,t) is 
the scalar charge density field arising from the presence of charged particles 
within the solution region, B(x,t) is the magnetic field within the region 
(which can be applied directly or induced by the flow of charged particles), 
and A(x,t) is the magnetic vector potential defined by (2.13f). 

Equation (2.13a), termed the Lorentz force equation, is the balance of 
particle inertia (on the left hand side) and particle accelerating force (on 


the right hand side). The time-dependent velocity component arises from 
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Newton’s second law of motion, namely 


dv 
F= = m— 2.14 


(where F is force and a the particle acceleration due to that force) and since 
dv /dt is the total derivative of the velocity field, we can write using the chain 
rule 

dv Ov 


pre Ca ear (2.15) 


so that (2.14) and (2.15) give rise to (2.13a). 

Equation (2.13b) represents conservation of charge and is a continuity 
equation indicating that the rate of change of charge density within a region 
is equal to the divergence of the current density from that region, where the 


current density J(x,t) is defined by 
J = pv. (2.16) 


Clearly if there is no net current density divergence from the region (no net 
charge flow from the region), then the charge density will remain constant 
within the region. 

Equation (2.13c) is Gauss’ Law and is the equivalent of (2.1) expressed 
in more than one dimension. It arises because sources of electrical charge 
can be mono-polar, and therefore give rise to a diverging electric field from 
a region enclosing them. 

Equation (2.13d) arises from the fact that magnetic monopoles do not ex- 
ist in nature, and hence there can be no net magnetic field divergence from 
a region. Therefore the magnetic field cannot consist of a scalar potential 
component, and can only be written in terms of the curl of a vector field as 


given by (2.13f). The link between the electric and magnetic fields is present 
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in the general electric field definition (2.13e). 


All the variables within the equations (2.13) above are time-dependent, since 
all fields and scalar quantities within the neutron tube acceleration region are 
time-dependent until an equilibrium state is reached. If the system is allowed 
to settle to such an equilibrium state (by applying constant drive conditions), 
it is considered steady, where no fields or scalar quantities change with time 


at any location within the region. 


2.3.1.2 Magnetic Field Influence 


At this point, the system of equations can be simplified by making the as- 
sumption that magnetic field effects are negligible. External magnetic fields 
are not applied to the neutron tube, and magnetic fields induced by the 
charged particle flow are small in comparison to the applied electric field. 
This is a known fact in the neutron generator community and is a conse- 
quence of the relatively low ion currents within the neutron tube; it is born 
out of extensive experimental experience and is justified to some extent, 


mathematically, in Appendix B. 
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2.3.1.3 Simplified Time-dependent and Steady-state Systems 


By removing magnetic field effects, the full system of equations, (2.13), can 


be restated as 


m i -V)v + a = gE, (2.17a) 
Op 

V - (pv) a (2.17b) 

V-Vo= af (2.17c) 
E0 


where in (2.17c), the electric field has been written as a conservative field 


(from (B.16) in Appendix B) with the definition 
E=-V¢4, (2.18) 


for the electric scalar potential ¢. This system of time-dependent equations 
can be further reduced, in the steady-state case, by immediately setting the 


time derivatives in (2.17a) and (2.17b) to zero, giving 


mv: V)v = -qV¢, (2.19a) 
V - (pv) =0, (2.19b) 
V-V¢= = (2.19¢) 

0 


In this case, the variables $(x), p(x), and v(x) are now the time-independent 
variables to be found. 

From a physical point of view, when in a steady-state, the ion current 
flowing across the tube acceleration region is constant, as are the consequent 
charge density, velocity and electric fields (provided that the applied voltage 
is constant). This leads to the immediate conclusion that the time derivatives 
within (2.17a) and (2.17b) are zero when the system is in this equilibrium 


state, hence (2.19), above. 


AT 


2.3.2 Solution Region and Boundary Conditions 


Depending upon whether a time-dependent or steady-state solution to the 
plasma boundary problem is required, solutions to the equations (2.17), or 
(2.19) are sought on the region R, shown in Figure 2.3, with the following 


boundary conditions applied 


o(x, t) = dg, t>0,x €S,US, (2.20a) 

o(x, t) = or, t>0,xEeS, (2.20b) 
V(x, t)-n, = 0, t>0,xES; (2.20c) 
V0(x,t)- ns =0, P0208 Ss (2.20d) 
P= pe b); t>0,x ES; (2.20e) 

p=, b> 0x € Soy (2.20f) 

p= Oe t=0,x ER, x €S,USz, (2.20g) 

V(x, t)- ns = 0, t> 0; x6 Ss (2.20h) 
v =0;(x)n, t>0,xES; (2.20i) 

y=, t=0,xER, x€S, (2.20}) 

v-n; =0, t>0,x €S; (2.20k) 


where {S;} are the boundary segments enclosing R (shown in Figure 2.3). 
The potentials @g and @r in (2.20a) and (2.20b) are the constant applied 
source and target voltages, respectively. The condition (2.20c) states that the 
electric field (given by (2.18)) normal to the plasma boundary (S}) is zero (ny 
being the vector normal to S;), and furthermore, since (2.20a) applies there, 
the electric field tangential to the boundary is also zero, or V¢(x,t)-s; = 0 


(where s; is the vector tangential to S,); this implies that Vé = 0 on S;. The 
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condition (2.20d) is applied to the line of symmetry, S;, and indicates that 
there is no change in potential across this boundary, whilst the additional 
condition (2.20c) on segment S, is required to determine the position of the 
plasma boundary, from which charged particles are emitted (|8]). 

Conditions (2.20e) and (2.20i) refer to the known charge density distribu- 
tion, p;, of ions at the emission boundary and the known velocity magnitude, 
v;, of the ions as they leave it (which is usually assumed to be constant in 
time). Ions emerge from the plasma boundary with a known kinetic energy, 
and it is assumed that they are emitted in a direction that is normal to the 
boundary surface, hence the presence of the unit normal vector nj in (2.20i). 
The condition (2.20f) states that there is no charge on the top of the plasma 
cup®, with (2.20g) and (2.20j) stating that the starting condition (for the 
time dependent problem) is a region without charge or velocity anywhere, 
except on those boundaries indicated. 

The final two conditions, (2.20h) and (2.20k), are applied to the line of 
symmetry, S;, and indicate that there is no change in charge density across, 
or that there is no velocity into the boundary, S;. The remaining segment 


does not have conditions applied. 


2.3.2.1 Overriding Problems 


The two overriding intrinsic problems with the model constructed here are 
firstly that from both the unsteady and steady systems of equations, it is clear 
the problem is non-linear in nature, and secondly that it is a moving (or free) 
boundary problem (we do not know the region upon which the problem is 


Tn reality this is a fixed metal boundary held at the specific electric potential dg. 
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to be solved in advance). Therefore, to understand the problem fully, the 
one-dimensional steady-state and time-dependent solutions are thoroughly 
studied, prior to examining a two-dimensional approach. As we shall see, 
the one-dimensional steady-state solution in planar geometry is analytically 
soluble, with the one-dimensional time-dependent case being at least par- 
tially, analytically soluble. Furthermore, studying numerical methods in one 
dimension allows us to check the validity and applicability of methods that 


might be applied to the two-dimensional case. 
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Chapter 3 


The One-dimensional 


Time-independent Solution 


The steady-state problem detailed in Chapter 2 can be restated in one di- 
mension. We now explore the steady-state one-dimensional case thoroughly 
with an view to gaining some insight into the one and two-dimensional time- 
dependent cases. As mentioned in §2.3.2, and as we shall see, the one- 
dimensional time-independent free boundary problem in planar geometry is 
analytically soluble. The solution to this problem is first derived, and com- 
pared with that given by a rapid numerical method before being restated in 
radial geometry, giving rise to a non-linear ordinary differential equation for 
which no closed solution appears to exist. The radial problem is then solved 
by modifying the numerical method used in the planar case (taking account 
of subtleties particular to the radial case), and solutions to this problem 
explored. 

In the following chapter, the full one-dimensional time-dependent solu- 


tion is examined and two different numerical approaches to this problem 
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described, along with numerical results. 


3.1 One-Dimensional Planar Time-independent 
Analytic Solution 


The equations, (2.19), can be written in an equivalent one-dimensional planar 


form as 


see ns al 
mo VT, (3.1a) 
d 
7g Pe) = 9 (3.1b) 

io p 

Ad = = (3.1c) 


where (3.1a) is simply the x component of the vector equation (2.19a), and 
m, q, and €9 are constants (as before). Here, what is now, the x direction 
was previously the vertical (y) direction. 

Intuitively, the continuity equation (3.1b) must be satisfied, since in one 
dimension, particles entering the solution region will always leave it (in more 
than one dimension this is not necessarily true due to the possible presence 
of magnetic fields causing vorticity within the charged particle beam). More 
importantly, (3.1b) implies that the quantity pv is constant at all points 


1 


within the one-dimensional region’. By integrating (3.1b) and setting the 


constant of integration to be Jo (the known initial current density) we have 


p(t) =, (3.2) 


'This must be so or pockets of local charge compression or rarefaction would un- 


physically form within the region. 
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or simply the one-dimensional equivalent definition of (2.16). 

The plasma boundary problem detailed in §2.3 is therefore reduced to the 
simultaneous solution for ¢(a) (the one-dimensional electric scalar potential 
or voltage) and u(x) (the scalar particle velocity) of (3.la) and (3.1c), and 
is considered on the domain x € [s, x7] (shown in Figure 3.1) subject to the 


conditions 


b= or, Top, (3.3a) 
a) as s; v= 8, (3.3b) 
db _ az 
or 0, Cas, (3.3) 
Y= Uj, L=s, (3.3d) 


where s is the location of the plasma free boundary. It should be pointed 
out here, that to accelerate a positively charged particle (such as a deuteron) 


away from the plasma boundary at s, then ¢s > @¢r. 


To solve the system, (3.1), we begin by writing (3.la) as 
d 2 
av (v(2) ) = aa 
which can be immediately integrated to give 
2 2q 
v(z)" = —— b(a) +e, (3.4) 


where c, is a constant of integration. When x = s, ¢ = ¢g and v(s) = 
uv; (the initial velocity of an emerging ion determined from the plasma ion 


temperature) from (3.3b) and (3.3d), so that the constant c; is given by 
2 
q= ve + ae 
m 
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Figure 3.1: The one-dimensional planar solution region showing the passage of the 
steady ion beam from the (settled) plasma boundary at x = s to the fixed target 
location at « = a7. Plasma expands from some location to the left of the plasma 
boundary and is accelerated away from the boundary by the electric field arising 
from the applied potential difference. The lines x = s and x = x7 correspond to 


the boundaries S; and S¢, respectively, in Figure 2.3. 


from (3.4). 


The scalar particle velocity v(x) can then be written 


v(x) =f! 5 — (2) +02, (3.5) 


since u(x) > 0 is consistent with particles travelling in the positive x direc- 
tion, which follows as ¢s > dr. Upon substituting (3.5) (along with (3.2)) 


into (3.1c) we have the single ordinary differential equation 

eb Jy (2q “3 

Ta = 2 (Ales ol) 402) (3.6) 
holding for s < x < xp and subject to the boundary conditions (3.3a), (3.3b), 
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and (3.3c). It is noted here that ¢ and a must be continuous on the interval 


[s, v7], whilst 9 must be continuous on the interval (s, x7). 


3.1.1 One-dimensional Planar Canonical Form 


In order to expose the underlying features of the one-dimensional problem, 
we now restate it in a simplified dimensionless canonical form, from which 
integration to give an analytic solution is relatively simple. 

Dividing the term in brackets on the right-hand side of (3.6) through by 
vu? gives 


i 
2 


Creo) +1) | (3.7) 


dx? EQui 


Oe Me Aa. 
mv? 


so that the dimensionless variable 


w(2) = Ls — (a) +1, 


> 0, (3.8) 


from (3.5), is the ratio of a particle’s kinetic energy at x to its initial kinetic 
energy. It readily follows from (3.7) that w satisfies 

a 

a = Bw-?, (3.9) 
(here it is assumed that w? = +,/w) where the constant, 

2qJo 

ae horses 

must have units of length, L. By defining y = (Gx and determining the 


necessary derivatives with respect to y, (3.9) can be further reduced to 


= p?, (3.10) 


which is the dimensionless canonical form of (3.6), where 


w(y) = w(y/B) = w(2). 


At this point, from (3.10), it is clear that w”(y) > 0, since 


from (3.8), and since v(z) must be positive for particles to travel in the 
correct direction. 
3.1.1.1 Dimensionless Boundary Conditions 


In order to solve the boundary value problem for (3.10), the original boundary 
conditions (3.3a), (3.3b) and (3.3c) must also be restated in a dimensionless 
form; this is done as follows. From the definition for w (or w) and the 


conditions (3.3a) and (3.3b), we can say 


= ko, (3.11) 


say, where yr = Grr. Clearly, ky > 0 (since dg > gr) and typically has a 
numerical value of the order 10° for the deuterium ions accelerated in this 


application. 


Defining 7 = (s consistently with the transformation from x to y, the re- 
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maining boundary conditions become 


Wyn = W(O(S)) 


=o (3.12) 
from (3.3b) and (3.8), and 
dw ; 
Ai = B*w'(s) 
Y ly=n 
-1_2q 
1 ! 
= -9 Ho's) 
=0, (3.13) 


from (3.3c) and (3.8). Similarly, as in (3.6), w and = must be continuous 


on the interval |y7, 7], whilst ab must be continuous on the interval (yr, 17). 
YT>7 ag 1 


3.1.2 Integration of the Planar Canonical Form 


The integration of (3.10) is easily performed by multiplying both sides of 
(3.10) by 252 giving 


di da dio 4 
dy dy? dy’ 
or 
d (dw? re 
— (=) =4—[2). 3.14 
dy (=) dy (w ) ( ) 


Integrating (3.14) we have 


Se 
om = 44/407 + 9, (3.15) 


where C2 is a constant of integration. The positive square root is taken due to 
the following argument. Since it has already been established that w”(y) > 0, 
then w’(y) is an increasing function of y for y > 7. Moreover, since w’(7) = 0, 


from (3.13), we conclude that w’(y) > 0 for y > 7. 
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Rearranging (3.15) at y = 7 and using the dimensionless conditions (3.12) 


and (3.13), we see at once that cp = —4 and therefore 
ey eee (3.16) 
dy 
If we now set 
w= wo? — 1, (3.17) 
then 
dw 
age (w? +1), (3.18) 
and since 
dw  dwdw 
dy dw dy’ 
then 
dw 
2(w?+1)—=1 3.19 
(w+ 1) (3.19) 


from (3.16) and (3.18). Integrating both sides of (3.19) gives the complete 
analytic solution of (3.10) subject to (3.11), (3.12) and (3.13) in terms of the 


variables w and y in the implicit form 


1 
where c3 is given by 
l ae 3 é 1 
c= 5ur— 3 (kb 1)" — (Bb 1)” a2) 


from (3.11) and (3.17). Finally, (3.20) can be written in terms of the original 
variables x and ¢ (from (3.8)) therefore giving the implicit one-dimensional 


steady-state analytic solution of (3.6) subject to (3.3a), (3.3b) and (3.3c) 


o(0) =5 ( 2 (¢5-0)+1-1) +( 24 (sa) 41-1] 
1 
> 


2 


me 


- (Gi 2 1), = (Ki = :)'] +2. (3.22) 


An expression for the location of the plasma boundary s for a given set 


of conditions can now be determined by setting « = s and consequently 


(8) = ds, giving 


_ 203 (3.23) 


from (3.21) and (3.22). 


3.1.3. Inversion of the One-dimensional Planar Solution 


The implicit solution (3.22) can be written in the explicit form ¢ = (x) by 
observing from (3.20) that, written in the transformed variables y = y(w), it 


is the canonical cubic polynomial. 


We can rewrite (3.20) as 


w? + 3wt+z=0, (3.24) 
where 
w= /(#? — 1), z=8(cs— 3). 


The discriminant D for (3.24), where 


DH]a108=977 


iy Vy, 


indicates that (3.24) has one real and a pair of complex conjugate roots, and 
since we are dealing with a physically real situation, we are concerned only 
with the one real root. To solve (3.24), we proceed in a manner analagous 


to Cardano’s method (circa 1545) and choose to detail the solution, instead 
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of using a symbolic solver. 


Equation (3.24) is similar to the identity 
4sinh® u + 3sinhu = sinh 3u, (3.25) 
and so writing w = 2sinhw in (3.24) we obtain 
2 (4sinh?® u -- 3 sinh w) +2z=0, 


or 


2sinh3u+ z= 0, 
from (3.25). Writing sinh 3u in exponential form we have 
eve "4+ 7=0, 
or by multiplying through by e®“, the quadratic 
(8)? + ze" —1=0. 


Since the exponential function takes only positive values, the appropriate 


solution to the quadratic is 
Qe" = —z + /(z? +4), (3.26) 


and by writing 
A=24+ f(z? +4), 


then 
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therefore 


N= e-8", (3.27) 


From (3.26) and (3.27) we conclude 
we (rztVE@t+O\E oa fztve +43 
eis 5 § Ser 5 : 
and that the inverted form of (3.24) is 
— (= v(2? oD) (= cs) 
I Oe ee epee et: calle 


This solution for w can be written in terms of the original variables x and 


2 
a 3(a-S)4 9(a- 2) +4 5 


which finally, upon squaring both sides and rearranging, gives the inverted 


form of (3.22) as 


— 13 (« — =) + 9 (« = =) + 4 7 (3.28) 
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where c3 is defined by (3.21). 


3.2 One-dimensional Radially Symmetric Time- 
independent Equation 


We now consider the plasma boundary problem on the domain shown in 
Figure 3.2. Here ions are accelerated radially outwards from the interface at 
r = s towards the target at r = ry, and their radial speed u(r, 0), the scalar 
potential distribution ¢(r, 0) and the boundary location r = s are determined 


by the boundary conditions 


b= or; r=rvr, (3.29a) 
o= os, Pas, (3.29b) 
dp 
aa 0, r=s (3.29c) 
U = Ui; res. (3.29d) 


As in the planar case, the potential ¢(r) decreases with increasing r, or 
és > or. Since the Dirichlet boundary value ¢r in (3.29a) is constant on the 
boundary at r = rr, and ¢g and v; in (3.29b) and (3.29d) are also constant 
on the boundary at r = s, then the scalar potential ¢(r, 0), ion speed v(r, @) 
(and consequent charge and current densities) are independent of angular 
position within the region and are therefore functions of radial location only 
(the plasma interface shape is also independent of angular location). The 
problem is therefore one-dimensional in r where r € [s,r7]. In the two- 
dimensional polar coordinate system, the gradient operator is given by 


ee hee. 
Va aah 
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Figure 3.2: One-dimensional radial solution region. The plasma is expanding 
from the centre of the two concentric circles (only half of the region is shown) and 
the plasma boundary (marked interface) forms some radial distance away from this 
central point at r = s. The target is located at r = rp. Since the ions are positively 


charged, dg > or. 


where p and q are the curvilinear unit vectors in the r and @ directions, the 


Laplacian operator is given by 


0 1 
vv ais (ts) + aoe 


and the divergence of a vector field F, say, is given by 


1d 10F, 
VR ae UFe) + ap 


where fF, and F, are the field components in the r and 6 directions, respec- 


tively. 
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One-dimensional radial equivalents to the system (2.19) are then found 


r Or Or r2 002 or Or Or 


to be 


OUp . 
=a P 
_ oo. 
=> 15, P> 
and 
10 1 Opv, 
Vv (NP Ue) TS 50 
+o (rpv,) 
Or ee 


since derivatives with respect to 0, and velocity components in the 6 (or q) 
direction are zero. To clarify, the one-dimensional radial equations for the 


plasma free boundary problem are 


oe = a (3.30a) 
ld 
oe (fpeye; (3.30b) 
ld /dd\ _ p 
aoe (3?) a, (3.30c) 


where the subscript p, denoting the radial direction, has been dropped and 
the partial derivatives replaced by ordinary derivatives. Here, v, ¢, and p are 
functions of r only. Manipulation of these equations can now proceed as in 


the planar case. 
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Integrating both sides of (3.30b) we have 


rpov = Ci, 
where c; is found to be 
Cy = SPsvj 
= sJo, 


ps being the initial charge density (at the free boundary), and Jp the initial 


current density (at the free boundary) from (2.16). We can then say 


J(r) = p(r)o(r) 


Ss 
= Jo- 3.31 
Jor, ( ) 


where the reciprocal dependence of J on r is expected for a radially expanding 
density. 

The speed v at r of a charged particle emitted from the free boundary at 
s can be determined from (3.30a) in a similar manner to the planar case (see 


(3.5) in §3.1 above) giving 


v(r) =f 4 (bs - (7) +02 (3.32) 


Using (3.31) and (3.32), we can now rewrite (3.30c), (3.30a), and (3.30b) as 


the single ordinary differential equation 


< (3?) ieee (= (bs - o(r)) + a) (3.33) 


éo \m 
where the unknown free boundary location s appears explicitly on the right- 
hand side. The solution of (3.33) subject to the boundary conditions (3.29a), 
(3.29b) and (3.29c) will result in a potential distribution as a function of ra- 


dial distance (that is particular to the specific boundary conditions applied). 
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3.2.1 One-dimensional Radial Canonical Form 


As in the planar case (§3.1.1), the ordinary differential equation (3.33) can be 
non-dimensionalised revealing a simple canonical form. Once again, dividing 
the term in brackets on the right-hand side of (3.33) through by v? gives 


f (3?) ee (“4 ($s — $(r)) + ne (3.34) 


2 
dr Eou; \ mu; 


so that the dimensionless variable 


= (3.35) 


can be defined, from (3.32). Clearly, since particles must accelerate away 
(in the sense of positive r) from the plasma boundary at r=s, then u(r) > 
v;, 7 > s indicating that w(r)2 Sahl Wee pee 


Derivatives of w(r) are then given by 


dw 2q do 
oe ph eee 3.36 
dr mv? dr ee) 
and 
du aq do 
dr2——— mv? dr?’ 
so that (3.34) can be written 
m2 ( dw _ dus sJo ot 
2 r— = = 
2q dr? dr EQu; b 
or 
d {fd 1 
ap (F) = sBw?, (3.37) 
where 
2qJo 
7 meu?’ 
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which has the dimensions of L~?. We now define the new dimensionless 


variables @ and R by 


and 
R= B2r, (3.38) 


respectively. Substituting these into (3.37) gives 


d dQ\ _ a 
a (x) =Q°2, (3.39) 


the dimensionless canonical form of (3.33), where Q(r) = Q(R/82) = Q(R). 


3.2.1.1 Radial Dimensionless Boundary Conditions 


As in the planar case, we must also restate the boundary conditions in a 


dimensionless form. With reference to (3.35) and (3.36) we have 


1 2 2 
=pris 8 ( “(bs - br) + i) 
mv; 
= 9-38 Fko, (3.40) 
say, from (3.29a) and (3.11), where Rr = Bare: 
Q(Rs) = Q(s) 
= B-3873 
from (3.29b), where R, = B28; and 
a] _ dQdu dol ar 
dR — dw dd dr|,dR 
R=Rs 
= 0, (3.41) 


from (3.29c) and (3.38). 
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3.3. Numerical Methods for the One-dimensional 
Plasma Boundary Problem 


There is no apparent analytic solution to (3.33) subject to the conditions 
(3.29a), (3.29b) and (3.29c). However, a number of numerical methods can 
be employed that utilise, for instance, moving mesh algorithms to iterate 
towards the solution of both (3.6) and (3.33) (an example is given in [34]). It 
is though, much more efficient to apply a relatively simple numerical method 
for initial value problems to these equations, than to employ a relatively 
complex moving mesh method. We now investigate the application of such 
a method to the one-dimensional, steady-state, planar problem, comparing 
numerical solutions to the analytic solution given in §3.1. This method is 


then used to study the one-dimensional, steady-state, radial problem. 


3.3.1 A Numerical Method for the One-dimensional Pla- 


nar Problem 


Since the one-dimensional planar equation (3.6) can be integrated analyti- 
cally giving the explicit solution (3.28), it would seem appropriate to develop 
a numerical method for this problem against which the analytic solution can 
be compared, before applying the method to the equivalent radial problem 
(3:33): 

A simple and efficient way of numerically solving (3.6) can be developed 
by splitting its canonical form, (3.10), into a coupled pair of first order ordi- 
nary differential equations, integrating them in non-dimensional space, and 


then returning the solution to physical space (where some physical inter- 
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pretation of the solutions can be extracted). This is most easily done by 


setting 
a (3.42) 
dy 
so that from (3.10) 
d 1 
a = 7}, (3.43) 


and integrating forwards with respect to y from the (source) free boundary 
position 7 = 3 2s (where s is selected in advance) with the initial conditions 
(3.12) and (3.13). Integration proceeds until (3.11) is satisfied, and since 
w'(y) > 0 for y > 7, the condition (3.11) can be satisfied only once over the 
solution region and this is at the right hand end of the region (integration 
proceeds with positive y from left to right). Therefore, the location of the 
right hand end (the target end) of the solution region is given by the value 
of y at the point where (3.11) is satisfied (where yr = Birr, say). 

The differential equation (3.10) is invariant under a translation of domain 
and this can be simply shown by shifting the dependent variable in (3.10) 
by a constant, a, say. Upon doing this, derivatives in (3.10) are unchanged 


yielding the similar one-dimensional equation 


1 


where m = y +a, and Wm = w(m). Here, the domain shift has introduced 
the second term on the left hand side of (3.45) 

With this in mind, it is clear that an initial choice of free boundary 
location 7 leads to a unique region size (yr — 7). Once this size is known 
(from the method described above), the fixed target location yr can be set, 
with the free boundary location at the left hand end of the region now varying 


with the applied boundary conditions. 
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3.3.1.1 Comparison Between Numerical and Analytic Solution 


To test the numerical method just described, solutions to the planar plasma 
boundary problem were found for several different sets of boundary con- 
ditions. The MATLAB function ODE45 was used to integrate (3.42) and 
(3.43) with constant initial particle energy (50eV) and current density (Jo = 
18800Am~*), typical of a neutron tube application, and with the target po- 
tential, @r, set at zero; several different source potentials (40, 60, 80 and 
100 kilovolts), that are typical of those seen in a neutron tube application, 
were applied. It was expected that the solution region would reduce in size 
as the applied potential difference was reduced (see Figures 1.2, 1.3 and 
1.4 and the associated text in §1.1.1), and this behaviour is clearly seen in 
Figure 3.3 which shows the numerical and the associated analytic solutions 
(calculated from (3.28)). The location of the free plasma boundary with an 
applied source potential of 100kV has been labelled as the nominal plasma 
boundary location to illustrate the reduction in region width with reducing 
source potential. The curve labelled Plasma boundary location as a function 
of source-target potential difference is calculated from (3.23). 

The ODE45 function was executed with a specified relative tolerance of 
1x10~°. This causes the solver, at a particular point, to repeatedly generate 
the solution with successively reduced integration step sizes. This step size 
reduction process is repeated until the difference between two consecutive ap- 
proximations of the solution, at that point, falls below the specified tolerance. 
In the planar case, this process can be checked by calculating the analytic 
potential solution, from (3.28), at the same integration locations as those 


generated in the numerical approximation; consequently the absolute RMS 
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difference (over the entire set of points) between the two solutions can be 
calculated. For the curves in Figure 3.3, calculated absolute RMS differences 
between numerical and analytic potential solutions range from ~ 7.5 x 107° 
to ~ 1.6 x 10~*. This gives confidence in the use of the ODE45 function 
so that it can be used in a numerical method for the subtly different radial 


problem, for which no analytic solution exists. 


3.3.2. A Numerical Method for the One-dimensional Ra- 


dial Problem 


There are some fundamental differences between the basic physical situations 
in the planar and radial cases that must be taken into consideration when 
calculating solutions to the radial problem described in §3.2. In the planar 
case, one of the initial specifications for the problem is the current density 
Jo at the ion source, with variations in this value giving rise to variations 
in the solution and solution domain size. The physically planar nature of 
the problem implies that the emitted ion current from the ion source is lin- 
early proportional to the emission ion current density, where the constant of 
proportionality is the emission surface area, and which in the planar case, is 
independent of domain location. This is manifested in the differential equa- 
tion (3.10), by it being invariant under a translation of domain; see (3.44) 
and associated text. Physically, this does not change the emitted ion current 
as already noted. 

The situation is different in the radial case where ions emitted with a 
given current density from a boundary close to the origin, give rise to a lower 


emitted ion current than those emitted from a boundary further away from 
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Figure 3.3: Numerical and analytic solutions to the one-dimensional planar plasma 
boundary problem at different applied source potentials (initial energy and current 
density is the same in each case). The accumulation of points towards the left of 
the calculated numerical solutions is caused by the error control procedure within 


the solver (described above). 


the origin. This is due to the emission boundary area, and hence emission 
current for a given current density, being radially dependent and is manifested 
in the differential equation (3.39), by it not being invariant under a shift of 


domain. This can be simply shown in the same manner as in the planar case 


re 


by shifting the dependent variable in (3.39) by a constant, a, say. Substituting 


the new variable, n= R +a, into (3.39), we deduce 


i (2) LO es (3.45) 


dn dn dn2 


where the domain shift has introduced the curvature term on the left hand 
side (not present in (3.39)), the effect on the solution of which, is clearly 
dependent upon the size of the domain shift; here Q, = Q(n). 

With this in mind, a particular domain location, size, and associated po- 
tential solution for this equation with a given emission ion current and set of 
initial conditions, is unique and cannot be translated (in order to fix the tar- 
get location) as in the planar case, hence a slightly different solution method 


must be employed and this is described as follows. 
To solve the radial problem described in §3.2, as in the planar case, we 
begin by splitting (3.39) in to a pair of first order equations by setting 
— =6, (3.46) 
so that expansion of the derivative on the left-hand side of (3.39) gives 
ee (@-4-<). (3.47) 


These can be used to calculate the single unique solution for a given set 
of initial conditions, by utilising the fact that Q(R) and Q’(R) are strictly 
increasing functions over the solution domain; this is shown in the following 


section. 
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3.3.2.1 Radial ‘Shooting’ Method 


As already mentioned, the situation in the radial case is subtly different to 
that of the planar case as a consequence of the plasma boundary location 
appearing explicitly on the right hand side of the original radial equation, 
(3.33). This is manifested in the dimensionless equation, (3.39), by it not 
being invariant under a shift of domain (unlike the planar case). This means 
that, in order to integrate forwards from the plasma boundary (as in the pla- 
nar case), the location of that boundary must be known in advance. Since 
this is not the case, a method for determining the boundary location, whilst 
simultaneously determining the solution that satisfies the boundary condi- 
tions, must be ascertained. 

In other words, selecting the plasma boundary location in advance and 
integrating forwards from it will not necessarily yield a solution whereby the 
boundary condition, (3.40), is satisfied at the required, fixed target location. 
Since integrating from the fixed target location is not feasible (there is no 
gradient information known in advance at the target), we must develop a 
method that iterates towards the correct plasma boundary location for a 
given set of conditions, whereby integration towards the target from the 
correct plasma boundary causes the boundary condition, (3.40), and the 
correct target location to coincide. This is reminiscent of a shooting method, 
commonly employed to solve ordinary differential equations (see [20], for 
example). 

To do this, we initially use a search method whereby we choose an initial 
plasma boundary location, R,, say, and and integrate towards the target, 


stopping when (3.40) is satisfied. If the integration stops before the required 
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target location is reached, the initial plasma boundary location is incre- 
mented by the small amount, AR,, and the integration repeated. This proce- 
dure stops at the 7“ initial plasma boundary location, Rs, = Rs,+(i-—1)ARg, 
when the chosen target location, Ry, falls within the range of calculated tar- 


get locations resulting from solutions starting at R 


s,_, and R,,. Once a 
plasma boundary location interval (containing the correct plasma boundary) 
that gives rise to a set of solutions, including that solution that gives rise 
to the correct target location for the given boundary conditions, is known, 
a bisection-type method is applied at the plasma boundary, with the calcu- 
lated target location being used at each step to inform the starting plasma 
boundary location; it being adjusted accordingly. 

This process is only possible if the dimensionless solution potential Q(R) 
monotonically increases (in a concave, upwards manner) from the known 
value at the plasma boundary to that at the target, or if both Q’(R) and 
Q"(R) are greater than zero over the domain R € [Rs, Rr]. To show this, we 


first make the transformation u = In R (so that u € [us, uo], with us = In R,), 


and Q(R) = Q(e") = P(u), so that (3.39) is reduced to 
P'(u) =e Po. (3.48) 


Now, since P(u) = Q(R), and 


>0, RE [R,, Rr] 
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(since 3, s > 0), we can conclude that 


‘ d dG 
P"(u) => RG (»%) 
aly (3.49) 


from (3.48). Additionally, by writing 


7 / e6P-2(¢) d¢ + P'(u,), (3.50) 


we can also conclude from the mean value theorem for integrals that P’(u) > 
0, since P’(u,) = 0, from (3.41); this, in turn, implies that Q’(R) > 0, R € 
[R;, Rr], from (3.50). 

It is now known that Q/(R) > 0 over the solution domain, from (3.50), 
but it is not immediately obvious that Q’(R) > 0, R €[Rs, Rr]. It is clear 
that the curve P(u) is concave upwards, since P’(u) > 0, from (3.49). Also, 
by superimposing the R axis over the u axis (but with a different scale), 
the graph of P(u) read against the R axis is simply Q(R); this must also 
be concave upwards as it is the same curve, merely being read against the 
different, strictly increasing scale, R = e". With this in mind, the curvature 
must be the same sign for both P(u) and Q(R), where for P(u), it is given 
by 


OP = Py 33 
(1+ (Pru))?)’ 
ou. u & [Us, Uo| 
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since P”(u) > 0. The equivalent expression in Q is then 
Q"(R) 
—— 
(1+ (@(R))) 


2Q = 


which has the same sign as gp, implying that Q”(R) > 0 (since Q'(R) > 0). 

Since it is the case that Q(R) curves upwards from R= R, in a concave 
manner, we expect that a change in initial plasma boundary location, for a 
given set of initial conditions, will give rise to a corresponding change (in 
the same direction) in calculated target location. This therefore enables the 


refinement process, described above, to take place. 


Radial Numerical Results 


As in the planar case, the MATLAB function ODE45 was used (again, with a 
specified relative tolerance of 1 x 107°) to solve the pair of equations, (3.46) 
and (3.47), thereby generating a potential solution for each initial plasma 
boundary location (as described above). The event trapping mechanism in 
ODE45 was used to terminate the integration when the value of the calculated 
(dimensionless) solution potential reached that given by (3.40). Dimension- 
less solutions were then transformed back to physical space for display. 

In each case, the initial search procedure starts the plasma boundary 
location guesses with the physical location sg = 0.001m, or dimensionless 
location R,, = 2.621. The search procedure successively increments from 
this starting location with an incremental step size of As = 5 x 10~+m, or 
AR, = 1.311. This continues until a plasma boundary interval is found, 
in which the corresponding, calculated target interval contains the required 


target location; the target was set at 0.02m in all cases. Once this interval 


1h 


is located, the bisection algorithm begins whereby the plasma boundary in- 
terval is successively halved (with the calculated target interval being used 
to inform the bisection method) until the calculated target location falls to 


within 10~'m (or 2.6 x 10~* dimensionless length units) of the required target 


location. 
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Figure 3.4: Numerical results from the radial shooting method. Curves show the 
variation in potential with distance away from the inner plasma free boundary (to 
the left) to the target (to the right). Solutions for a number of differing sets of 


initial conditions are shown. 
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Figure 3.4 shows results from two current (1 Amp and 10 Amps) régimes 
calculated by the radial shooting method just described. As is expected, cal- 
culated plasma boundaries within the device carrying the smaller 1 Amp ion 
current are located nearer the origin than those within the device carrying 
the larger 10 Amp ion current, for a given external electric field. This is due 
to the larger ion flux in the 10 Amp case causing the boundary to bulge into 
the accelerating region (and conversely for the smaller 1 Amp case) as de- 
scribed previously in 81.1.1. The effect of radial location on the solution can 
be clearly seen in the 1 Amp case, particularly for solutions with a (relatively) 
high applied potential difference. The combination of a strong electric field 
and relatively low ion flux causes the plasma boundary to be located near to 
the origin, where the circular boundary begins to resemble a point source. 
This is in contrast to the planar case and is manifested in the presence of s 
on the right-hand side of (3.33), in addition to the more complex derivative 
term on the left-hand side of (3.33). For certain sets of boundary conditions, 
these terms cause the solution curvature to change sign (from negative to 
positive) some distance away from the plasma boundary. 

In the planar case (3.6), the electric field gradient (and hence curvature) is 
always less than or equal to zero, irrespective of the location of the free plasma 
boundary. However, in the radial case (3.33), the electric field gradient is 
given by 

@d 1fsJy ( 2¢ “2d 
aa e (= (4 (és - 0(7) +1) +2), 
where, since a <0, r € (s,rr] (from (3.35) and (3.36)), we conclude that 


the solution curvature changes sign where 
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(3.51) 


Referring to the 100KV, 1 Amp case in Figure 3.4, the solution gradient 


is roughly constant, taking a value of ~ —6.25 x 10° Vm~'. In combination 


with the constants 


and, 


2q il 
mue Ip 
=0020 J? 
sJo Io 
Egui ~ TSE QU; 


d 


~ 1.29 x 10’ C?(Jm)~’, 


where I, is the initial particle energy in electron-volts (Iz = 50 eV in this 


case), the current density Jo is calculated as the density over the circular area 


of the emission boundary (Jo = Jo/(7s”) ~ 210 Am~”), and s is the plasma 


boundary location (~ 0.04 m in this case); (3.51) indicates that the solution 


curvature will change sign where 


os — o(r) < 21500V. 


(3.52) 


A simple visual inspection of the 100KV, 1 Amp trace in Figure 3.4 indicates 


that the solution curvature changes from sign when the solution potential is 


at ~ 80000V (or where dg — ¢(r) ~ 20000V), corresponding to (3.52) to a 


first approximation. 
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Chapter 4 


The One-dimensional 


Time-dependent Solution. 


4.1 Introduction to the One-dimensional Time- 
dependent Problem 


As mentioned in 81.1.1, two distinct timescales exist during the operation 
of the neutron tube, and the previous chapter described (in one dimension) 
the stable steady state to which the system settles over the relatively long 
timescale. In this chapter we examine (in one dimension) the relatively short 
timescales involved with the passage of the charged particle beam across the 
tube accelerating gap, and the consequent settling of the plasma-vacuum 
boundary. 

Prior to the creation of the deuterium plasma, a state exists within the 
operation of the neutron tube, where no charged particle flow takes place. In 


this state a stable electric field is generated by the application of the main 
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acceleration voltage to the tube electrodes. As the deuterium plasma is re- 
leased, ions begin to accelerate across the tube accelerating gap separating 
the solution region into two states. In the part of the region ahead of the 
advancing charge front there is clearly no charge (i.e. the charge density and 
velocity fields are zero), but behind it the region is charge infused, where by 
charge infused, we mean that this part of the region contains charged par- 
ticles, with the charge density and velocity fields being non-zero. Since two 
states exist within the solution region as the charge front is accelerated across 
it, the solution of the time-dependent (short timescale) problem can be anal- 
ysed differently in these two areas. The charge free region ahead of the charge 
front can be analysed analytically, as we shall see, with it being effectively 
separated from the charge infused region by a characteristic curve of the 
problem. This “separation” characteristic is defined by the time-dependent 
location of the charge front itself. A partial time-dependent analytic solution 
can be also obtained for the charge infused region behind the charge front, 
but to obtain a full time-dependent solution both behind and ahead of the 
charge front, we appeal to a numerical method. 

Two numerical approaches are proposed for determining the solution of 
the time-dependent problem, where the first numerical method, that uses 
a standard difference approach, required, in hindsight, a parameter to be 
introduced to overcome stability issues that are inherent to the method and 
that lead to its failure. A new, second numerical approach, was devised 
to overcome the difficulties experienced with the first method, and this is 


detailed in Section §4.3.3. 
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4.1.1 Fields Pertinent to the Problem 


0=Os 
~~ """“Advancing wavefront. 


o 
S -S- 
a I 
S ae 
: Charge infused Charge free 
region region 
x=s(t) x= 


Figure 4.1: One-dimensional planar solution region showing (in the short 
timescale) the rapidly advancing charge wavefront. Plasma expands from some 
location to the left of the time-dependent plasma boundary (at «2 = s(t)) and is 
accelerated away from the boundary by the electric field arising from the applied 
potential difference. The region to the left of the wavefront has non-zero charge 
density, whilst that to the right has zero charge density. The target is fixed and 


located at x = x7, as before. 


At this point it is valuable to summarise the various vector and scalar 
fields pertinent to the time-dependent problem before restating the solution 
region and time-dependent system of equations in a one-dimensional form. 
We assumed in §2.3.1 that magnetic field effects can be ignored and so the 


relevant variables are as follows. 
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The charge density scalar field is represented by p(x,t), and represents 
the charge per unit volume within an enclosed region. It is related to both 
the current density and velocity vector fields via the definition (2.16), where 
the current density J(x,t) is the time-dependent rate of flow of charge per 
unit cross sectional area at x and the velocity v(x, t) is the time-dependent 
rate and direction of movement of charge at x. Finally, the electric field 
E(x,t), defined by (2.18), is the time-dependent force acting on a unit of 
charge at x. 

In one dimension, all fields are spatially dependent upon only one vari- 
able, x say, and Figure (4.1) shows the one-dimensional solution region with 
the charge front (separating the two charge infused and charge free regions) 
advancing in a positive x direction, from the plasma boundary on the left at 


s(t), to the target on the right at xr. 


4.1.2 The One-dimensional Time-dependent System 


The time-dependent system, (2.17), in now restated in one spatial dimension 


(the spatial independent variable being x) as 


Ov Ov’ q 
O(pv) _ Op 
A Oe (4.1b) 
ro p 
ee) = oe (4.1c) 


in addition to the one-dimensional conservative electric field definition, 


et, (4.1d) 
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The required time-dependent solution is generated by solving these equations 


on the region s(t) <a < xp and t > 0, in conjunction with the conditions 


SG hi belie t>0 (4.2a) 
Md be oie t>0 (4.2b) 
2ets(.0) =0, t>0 (4.2c) 
COR ones t>0 (4.24) 
0: PC (4.2e) 
p(s(t),t) = ps(t), t20 (4.2f) 
Asse $e Gee (4.28) 


Here s(t) is the, as yet, undetermined time-dependent plasma boundary lo- 
cation with sp = s(0) = 0 (so is set to zero in this chapter, but need not be 
so); @r, ds, and v; are known constants with ¢s > @r; the function p,(t) is 
a known time-dependent function, which could take any physically realistic! 
form, but is known from experiment to be a function ramping from zero to a 
constant value over a relatively short period of time (being constant for the 
majority of the neutron tube operation time). 

The advancing charge density wavefront is represented by the coupled 
hyperbolic equations (4.la) and (4.1b), whilst the electric field in both charge 
infused and charge free regions is coupled to these equations by the elliptic 
equation (4.1c). The electric field, in the absence of time varying magnetic 
fields is defined as (4.1d), the gradient of some scalar potential ¢, as in (2.18). 


'By physically realistic, we mean a function that can be driven by an experimental 


procedure; one that is not discontinuous. 


85 


4.2 Analytic Solution 


By analysing the time-dependent system, (4.1), inroads can be made in to an 
analytic time-dependent solution. Initially, we examine the solution ahead 
of the advancing charge front, and follow with further analysis that gives rise 
to the analytic form of one variable, within the full time-dependent prob- 
lem, upon characteristic curves pertinent to the hyperbolic equations, (4.1a) 
and (4.1b). This analysis reveals, by manipulating the original system of 
equations, a second order Riccati equation that can be linearised and solved 
yielding the explicit form of the variable 0v/Ox upon characteristic curves 


spanning the solution region behind the charge front. 


4.2.1. Analytic Solution Ahead of the Advancing Charge 


Front and Associated Separating Characteristic 


As previously described, when charged particles are introduced into the ac- 
celerating region, they propagate across it in a wavelike manner, effectively 
partitioning the region into two areas that can be considered separately. 
The part of the solution region ahead of the propagating charge front is 
charge free, whilst that behind it is charge infused. Within the charge infused 
part of the solution region, the electric field is inherently dependent upon the 
charge density field as a consequence of Gauss’ Law, (4.1c), and is therefore 
time-dependent, since the charge density is time-dependent as a consequence 
of the conservation law (4.1b). However, ahead of the advancing wavefront 
the electric field is considered to be constant as a consequence of the absence 


of charged particles”. The electric field and consequently the charge density 


Strictly, as a consequence of the infinite range of the electromagnetic field and of its 
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field are coupled to the velocity field via the Lorentz force equation, (4.1a). 

At t = 0, no charge exists within the acceleration region as no particles 
have been introduced at this time; this is formalised with condition (4.2g). 
Clearly, due to the absence of charge across the acceleration region at t = 0, 


(4.1b) is automatically satisfied. Furthermore, with p = 0, (4.1c) becomes 
ro 
Or? 

which upon integrating and applying the conditions (4.2a) and (4.2b) gives 


the expression 
bs — Or 
IT — 80 


o(e,0) = ( 


for the electric scalar potential within the acceleration region at t = 0. Also, 


) r= 2) + or (4.3) 


from the definition (4.1d) and (4.3), the electric field E(x,0) must be the 


constant 


=e (Sa 2) (4.4) 


LT — So 
An expression for the velocity of particles emitted from the t = 0 boundary 
can now be obtained by firstly writing (4.la) as the ordinary differential 
equation 


Ni 4. 
Ho me” (4.5) 


where the time, t, can be measured along characteristic curves in the (z,t) 


plane. 

speed of propagation, c (the speed of light), it is conceivable that the low density presence 
of charge at the advancing wavefront could influence the electric field ahead of it, with 
this influence decaying with distance from the charge front. Such an influence is believed 
to be very small and is therefore assumed to be zero, with the electric field ahead of the 


advancing front assumed to be constant. 
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Integrating (4.5) we have 
[qunt feu, 
dt m 


v= 2 Bet + Ca, (4.6) 
m 


or 


from (4.4), giving rise to the following expressions for the velocity field along 


curves originating at t = 0. Applying the conditions (4.2d) and (4.2e) we 


have 
v= Iey + U; (So 0) (4.7) 
m 2 2 ? 
v= = Bt, (2,0), (4.8) 


where (so, 0) represents the origin of the curve at (x = s9,t = 0), and (2p, 0) 
the origin of the curves at (v = x»,t = 0), and where x, € (so, 27]. Ex- 
pressions for the family of characteristic curves can now be determined by 


writing (4.6) as 
dx q 
a a 
dt m° aoe 


which, upon integration, gives rise to the curves 


5 Bot? + uit + So, (So, 0), (4.9a) 
q 
‘i 5, ot” + Zp, (35:0); (4.9b) 


from (4.7) and (4.8); here, as before, (so,0) and (x,,0) represent the origins 
of the curves. 

Equations (4.7) and (4.9a) are simply the expected speed and position at 
time t of a single particle (within the acceleration gap) with a charge-mass 


ratio of (g/m) under a uniform acceleration of magnitude + Eo, emitted from 
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the initial plasma boundary (located at s9). Additionally, whilst (4.9b) and 
(4.8) represent the location and gradient (respectively) of characteristics in 
the charge free part of the solution region, they also represent the expected 
speeds and positions of single particles, if they were emitted from differing 
locations (x,) between the initial plasma boundary (so) and the fixed tar- 
get (xr). The characteristic curve (4.9a) also separates the charge free and 
charge infused parts of the solution region, and represents the location of the 
advancing particle wavefront. Figure 4.2 shows the curve (4.9a) and family 
of traces (4.9b) for single deuterons within an acceleration region bounded 
by a fixed target located 0.01m from the initial plasma boundary, and with 
an applied potential difference of 120kV. The part of the solution region 
enclosed by the separation characteristic, the lower t = 0 boundary, and the 
right hand fixed target is charge free, with the electric field being constant 


(Eo) everywhere within it. 


4.2.2 Analytic Solution Behind the Advancing Charge 


Front 


A further analysis of the time dependent system, (4.1), has been developed 
and is now described. This new analytic approach to the problem, once 
again, considers suitable characteristic curves giving rise to a partial analytic 
solution to the full time-dependent problem; it yields the closed form of the 
variable Ov/Ox on characteristic curves behind the advancing charge front. 


The analysis is as follows. 
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Figure 4.2: Actual characteristic curves spanning the charge free part of the so- 
lution region enclosed by the separation characteristic, the lower t = 0 boundary 
and the right hand fixed target. An impression of the moving plasma boundary is 


shown (marked s(t)). 


By taking the derivative of (4.1a) with respect to x, we have 


dv\? 8 (du da (dv q OE 


Once again, introducing a family of characteristics, where in this case it is 
useful to introduce the coordinate o, so that x = x(o), t = t(a), with o being 


measured along the characteristics and being defined by 


dx 
a 4.11 
dt 
—-=1 4.12 
da ? ( ) 
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allows (4.10) to be written on the characteristics as 


d (Ov\_ q dv \? 
de (=) See (a) ry 


from (4.1c); here it is assumed that 


a (20) _ a (a 
Ox \ Ot) Ot \Ox)~ 


Defining the function 


and writing (4.13) in terms of a gives 


da q 2 
— = — p- a’. 4.14 
do Wee = ( ) 
Now expanding (4.1b) we have 
Op Op Ov 
fal rect igen 4.1 
"an Ot P On’ ee 


and when similarly defining p(o) = p(x(c),t(c)) on the characteristics we 


have 
dp 
—=- 4.16 
qo OM (4.16) 
for (4.15). Taking the derivative of (4.14) with respect to o gives 
d (da 2\_ @ dp 
do Gaz ) eae 
MEQ , 
from (4.16), or, 
1@a da q 
——~ + 2 + —p=0. 4.17 
a do? da o még’ ( ) 


Substituting (4.14) into (4.17) gives the second order non-linear ordinary 


differential equation for a(c) 
d 
ee Be va a (4.18) 
o 
on the characteristics. 


4.2.2.1 Linearisation of the Riccati equation 


The non-linear ordinary differential equation (4.18) is a Riccati equation of 


order two, and can be linearised by making the substitution ([25], [35], [43]) 


1 dz(c) 
= —~ 4.1 
a(a) = zo), (4.19) 

whereby 

da 1d@z 1 (dz\? 

ee ee | 4.2 

ag Zao? ~ 2? (=) eee 
and 


2 3 2 3 
dia _ldz 3dedz | 2 (=) . (4.21) 


de® zdo® dodo? 2 \do 
Substituting (4.19), (4.20) and (4.21) into (4.18) reduces it to the particularly 


simple third order equation 
which can be immediately integrated giving 


z(o) = Ao? + Bo + C, 


3Tt also known as a modified Emden-type equation, or modified Painlevé-Ince equation. 


v2 


where A, B, and C are integration constants. The variable a(a) can now be 


recovered giving 


Ov 
Ox |, 
2Ao0+ B 


Ao? + Bao+C 

20+ By 
= 4,22 
o2+ Bio +C,’ ( ) 


from (4.19), where B; = B/A and C, = C/A. From (4.12), we see that 


a(o) = 


go =t-—to, where top is a parameter representing the origin of that particular 


characteristic on the s(t) boundary, so that a(o) can be written 


a(t) = a(a(t)) 
Ov 
ae |, 


2t + By 
Ge eet 4.23 
2+ Bot + Cy’ oe) 


where By = By, — 2to, and Cy = C; + t2 — Byto, from (4.22). 
4.2.2.2. Determination of 0v/Ox at the plasma boundary. 


Since v = v; at the plasma boundary (from (4.2d)), then its derivative (fol- 


lowing the boundary) is zero there, or 


yy OU. 100 
Additionally, we have from (4.1a) 
Ov Ov 
fuse + a Bors =a 0, (4.25) 


since E(s(t),t) = 0 from (4.2c), and again noting v(s(t),t) = v;. Subtracting 
(4.24) from (4.25), 
Ov 
s=s'(t))=— =a) 
le-rog] =o 
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and since vu; > s‘(t) (or particles would never be able to leave the plasma 


boundary), we conclude 


x 
a =, (4.26) 
i: x=s(t) 


By defining the time t, to be the time when the plasma boundary is at the 
location s(t,), then (s(t;),t;) is the origin of a particular characteristic on 
the plasma boundary. It is possible to deduce that characteristics originating 
at the plasma boundary monotonically increase without crossing each other, 
so that t = t, once, and only once, on a particular characteristic. With this 
in mind, we can apply the boundary condition (4.26) to (4.23), thereby de- 


termining the constant By in terms of t,. 


It is also possible that the constant C2 can be determined, although the 


resulting equation for Ov/Qz is likely to be very implicit. 


4.3. Numerical Schemes for the Determination 
of the Time-dependent Solution 


Progression to a full analytic time-dependent solution appears difficult, and 
so we now appeal to numerical algorithms for the determination of the full 
time-dependent solution of the system, describing the two numerical ap- 
proaches mentioned above. 

Sets of results have been chosen to highlight features of the two methods, 


which were both implemented in MATLAB. 


The physical conditions generally applied to the problem are listed as fol- 
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lows: 
e General Physical Conditions for the Problem 


— Target potential @r = 100 kV; 
— Source emission potential ds = 0V; 


— Particle emission velocity v; = 69181.8 ms~! (deuterons emitted 


with an energy of 50eV); 
— Source current density Jo = 1.88 x 10* Am7?; 


— Initial region size L = 0.01 m. 


4.3.1 Reduction of the Time-dependent System to Di- 


mensionless Form 
Before the numerical methods are developed, the system, (4.1), and associ- 


ated boundary conditions, (4.2), are re-written in dimensionless form. To do 


this, we can define the dimensionless variables 


d= -. (4.27a) 

p= - (4.27b) 

es U 

v= ie (4.27c) 

z= a (4.274) 

E= = (4.27e) 
and, t= a (4.27f) 
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where @, = ds — ér is the potential difference across the region, p;, is the 
maximum charge density (naturally located at the moving emission bound- 
ary), Um is the maximum particle speed (located at the target boundary), 
L = xr — 8p is the initial region size, Ey = ¢,,/L is the constant initial 
electric field magnitude, and t,, the estimated time taken for a particle to 
traverse the region. For the purposes of developing the dimensionless set 
of equations, it is convenient to assume that particles are emitted from the 
emission plasma boundary with no initial energy, even though this is not 
true in reality. With this in mind, a particle reaching the target, having 
undergone linear acceleration across the region from an initial speed of zero, 


will have a final, maximum speed 


Um = \/ 240m (4.28) 
m 


where q is the particle charge and m its mass, as usual. At this speed, a 


particle will then take the time 
tn = — (4.29) 


to cross the initial region, and we shall use t,, as the time scale. 
Using these variables, the system, (4.1), can now be written in dimen- 
sionless form. To do this, initially, the variables (4.27a), (4.27d), and (4.27e) 


are substituted into the electric scalar potential definition, (4.1d), giving 


Dean cy (4.30) 
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Also, substituting the variables (4.27b), (4.27d) and (4.27e) into (4.1c) gives 


OE _ Epa. 
0% 7 Pest 
__&6 
Ox?’ 
from (4.30), or 
ad 
ap OP: (4.31) 
where the dimensionless parameter, 
ogre 
= 7 4.32 
7 PmEo ( ) 


has a typical value of O(10). 
Additionally, substituting (4.27c), (4.27d), (4.27e), and (4.27f) into (4. 1a), 


we have 7 
Oo 6 _ Ov 


a 2 dR’ 
and finally, substituting (4.27b), (4.27c), (4.27d), and (4.27f) into (4.1b) we 


(4.33) 


have 


06 Op 
we & + ise , (4.34) 


od e 

as = —a1/p, (4.35a) 
aap ie, (4.35b) 
Op 00  _Op 
a Ge eb: ch : (4.35c) 


which are to be solved on the region §(t) < £ < Zr, where 


soe SUE) 
2 x 
and, ¢7= A 


The equations, (4.35), are to be solved in conjunction with the dimensionless 


boundary conditions 


dr, t) = 2 
= ¢r=const, £>0 (4.36a) 
5(a(),8) = = 
= $5 =const, t>0 (4.36b) 
2 a ==0, t>0 (4.36c) 
1(8@).8) = 
=%;=const, t>0 (4.36d) 
3(#,0) = 0, # € (8, &r] (4.36e) 
if), f) = 80) 
= p,(t), t>0 (4.36f) 
p(#,0) = 0 & € (8, &r], (4.362) 


where S9 = S(0). It is noted that if dr = 0 (which is often the case in 


experiment), then dr = 0 and gs = 1. 
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4.3.2 Numerical Method 1 for the Determination of the 


Time-dependent Solution 


The three dimensionless, coupled equations, (4.35), can be used in a time 
stepping algorithm to attempt to determine the time-dependent solution to 
the plasma boundary problem in one dimension. The proposed algorithm is 


broadly as follows: 


1. The initial solution region is set with the initial plasma boundary lo- 


cation being defined as 5(0) = 50; 


2. Set p(Z,0) = 0(@,0) = 0, where Z € (5, Zr]; 


3. Set p(5,t) = p(t) and 0(5,t) = 0;, from (4.36d) and (4.36f), where 
fs(t) is the dimensionless time-dependent function representing the 
time-dependent variation in current density at the plasma boundary 
during tube operation, and 0; the initial dimensionless particle velocity 


(a constant in this case, although not necessarily so); 


4. Calculate the potential field 4(2, t*) from (4.35a), using the conditions 
(4.36a) and (4.36b), and the most recent value of the charge density 
field (t* refers to the time at the k'” time step with ¢° being the initial 


time step; this is detailed later); 


5. Calculate the electric field E(é, t*) from the potential field b(é, t*), 
using (4.30); 


6. Calculate the updated velocity field, 0(z,t*+!), from (4.35b), using 
E(é, t*), and 6(&, #*); 
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7. Calculate the updated charge density field, p(z,t*+!), from (4.35c), 


using A(Z,t*) and 6(z, t**+) from Step 5; 


8. Investigate the electric field magnitude E (3,t**") at the present loca- 
tion of the emission boundary. If this is not zero, by calculating the 
electric field gradient and curvature at this point, predict the location 
of the boundary where it is likely to be zero and adjust the boundary 


location accordingly; 


9. Return to step (3) and repeat until the solutions in successive time 


steps are identical to within a small assigned parameter. 


4.3.2.1. Mapping from the Moving Physical Domain to a Fixed 


Logical Domain 


As a consequence of the moving plasma boundary, S(t), the solution region is 
time-dependent. If the system of equations, (4.35), is to be discretised using 
differences (as is proposed here), then nodal spacing in the spatial domain will 
be time dependent, introducing a variation in accuracy at each time-step. In 
order to avoid this, we can map the differential equations, (4.35), to a fixed, 
logical region with a logical spatial independent variable € € [0, 1], and logical 
time variable 7 € [0,00). By doing this, since the physical domain is time- 
dependent, the mapped differential equations will change in form with time, 
such that a different system is solved in logical space at each point in time. 
However, since the logical region is fixed, equally spaced differences can be 
used to solve the different differential equations arising as a consequence of 
the mapping at each time step, with the physical solution being reconstructed 


subsequently. 
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To perform the mapping, the dimensionless physical independent. vari- 
ables are now written % = Z(€,7), and t = ¢(€,7), and in order to introduce 
a method of varying the physical spatial variable x with the logical variable 
in a non-linear, time-dependent way, we introduce the mapping function 
M(€,7). Here, we choose 
) = a 

0g 


K5(T) 


SS SE (4.37) 


Vl + (fo? 


where the function f(E, tT) = f(&(E,7), t(€,7)) is a chosen required physical 


M(é,7 


solution variable mapped to the logical region (the subscript variable de- 
noting a partial derivative with respect to it and two subscripted variables 
denoting the second partial derivative etc.); 4 is a chosen parameter. The 
mapping function, M(€,7), is similar to monitor functions that are com- 
monly employed in moving mesh methods (see [33], for example); a detailed 
analysis of their use can be found in [42]. By using a mapping function of the 
form (4.37), the spatial nodal density can be made to adapt to the solution 
at the previous time step and so the use of (4.37) provides a mechanism for 
adaptively controlling the spatial nodal density (also used in [30]). 

The scaling parameter, &,(7), is a time-dependent normalisation param- 
eter used to define the physical region size (since the physical region size is 
time-dependent); it must be extracted as part of the solution, to determine 
the physical nodal density at a particular time. Clearly, where i (€,7) varies 
rapidly with €, then depending upon the magnitude of the parameter pu, the 
mapping function M(£,7) (and hence the rate of change of physical vari- 


able with logical variable) is small (thus in a discrete scheme, nodes will be 
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concentrated in this area). Where f(€,7) varies very slowly (or not at: all) 
with €, the mapping function M(£,7) is approximately constant in €, and 
is equal to the scaling parameter «,(7). The parameter js can be chosen to 
offer further control over the rate of change of physical variable with logical 
variable. 

Similarly, in addition to (4.37), and since we expect the plasma boundary 
position and hence the solution region to be rapidly changing, we control the 
rate of change physical time ¢ with logical time 7 by introducing the time 
step controlling function G(r) given by 

a = G(r) 


Ct 


= = (4.38) 


where c and vy are chosen constants, and g(7T) is a suitably chosen time- 
dependent function, the derivative of which, g,, is significant during times 
where the region size is changing rapidly. When the region size and conse- 
quently the function g(7) varies rapidly with time, the function G(r), and 
hence the rate of change of physical time with logical time, is small. Dur- 
ing periods of time where the region size is not changing rapidly, the rate 
of change of physical time with logical time will reduce to the constant c;. 
Again, the parameter v can be chosen to offer further control over the rate 
of change of physical variable with logical variable. 

We can now re-write the system of equations (4.la) to (4.1d) in terms of 
the logical variables € and 7, by firstly noting that we choose the physical 
time variable t to be independent of the logical spatial variable, €, or that 
¢ = t(r) only. With this in mind, the derivatives within the system are trans- 


formed as follows. 
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Using the chain rule 


oad, aa 
OE O06 0% OE Ot 
0 
= M— 4, 
from (4.37), noting that 0#/0€ = 0, or 
0 10 
Also 
OF PO 0 
Or Mdé Ot’ 
from (4.38) and (4.39), or 
Oi. WP Oe «710 O 
a= Gar” rae): oa 
Furthermore 
oP 1 1 0M O 
am MPOE MP OE OE" aa: 


from (4.40). The system of equations can now be restated in the logical do- 


main as follows. 


Using the above transformations, the non-dimensional definition (4.30) can 


be written as 


E(#(E,7), t(€,7)) = EE, 7) 


= (4.43) 
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from (4.40), whilst (4.35a) can be rewritten 
oe 106 1 AMAd 


O72 M2 0&2 M3 OE OE 


= —a1p, 


from (4.42), where (€,7) = (&(€,7), t(€,7)). 


(4.44) 


Furthermore, using both (4.40) and (4.41), the time-dependent equations, 


(4.35b) and (4.35c) can be similarly transformed giving 


ue i ee (eee 
de Ot M G/ @€ © Gar 
_E 
=> 
or 
a E 17¢ a 
aay ea a (s Fat 
and 


or, on expanding and rearranging 


G{.06 7, £\Odp) a 
Pag * (? =) 9¢} = Or’ 


where 0(#(€,7), #(€,7)) = O(E,7). 


(4.45) 


(4.46) 


(4.47) 


(4.48) 


Now that the system of equations has been rewritten in terms of the 


logical variables € and 7, it can be solved on the fixed logical domain by 


noting that the conditions at the boundaries of the non-dimensional physical 


domain, (4.36a) to (4.36g), also apply at the boundaries of the logical domain. 


This is because, for example, 


(4.49) 


and thus the potential condition (4.36a) on the target at Zp is then 


It is noted that the derivative condition (4.36c), applied to the plasma bound- 
ary in the physical region, is not explicitly applied in the logical region. This 
condition is utilised in Step 8 of the algorithm above, where it is used to 
determine the time-dependent location of the plasma boundary in the phys- 
ical region. Since the plasma boundary location is adjusted in the physical 
region at each point in time, the physical region size and hence the forms of 
(4.43), (4.44), (4.47) and (4.48) are time-dependent, due to «,(7) in (4.37) 


(as previously mentioned). 


4.3.2.2 Numerical Implementation of the Time-stepping Algorithm 


Using Differences 


In order to solve the system (4.43) to (4.48) on the logical solution domain us- 
ing the algorithm above, we initially split the domain into n discrete, equally 
spaced subregions, giving rise to the n + 1 nodes {€;} € [0,1] separating 
the subregions. Each logical subregion has the size A€ = 1/n, so that the 
j node has the logical location €; = (j — 1)A€, where j € [1,n + 1] and 
j € Z. Additionally, within the problem, the time 7 must fall in the range 
T € (0, co), with 7 = 0 indicating the point at which ions are first accelerated 
away from the emission plasma boundary. Representing time in a similarly 
discrete way, the k** time step (with k € [1,oo) and k € Z) can then be 


written 7", where {r*} € [0, 00). By choosing a fixed time step size Ar, the 
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k** time step is then at the time r* = (k — 1)Ar. Using this notation, the 
value pr is the value of the charge density, say, within the logical region at 
the location (7 — 1)A€é and time (k — 1)Ar. 

To determine solutions to (4.43), (4.44), (4.46), and (4.48) on the logical 
domain, we must first determine an expression for the time-dependent scaling 


parameter «,(7). This can be done by integrating (4.37) to give 


ip = 3() 4 : M(é,r) dé 


ao fa 
Vvit+pe co 


for the n + 1 nodal points {€;}, where 8(¢) is the time-dependent plasma 


= 5(t) +6 (4.50) 


boundary location; #7 is the fixed target location. Rearranging (4.50) gives 


the required expression for k,(T) as 
= 


41 
Ks(T) = (Zr — s(t a aac 


Physical node locations (and hence the solution in the physical domain) can 


(4.51) 


then be recovered from (4.50), with the p“” nodal location being given by 


€j+1 dé 
pay) dae Vie 


Similarly, discrete time points in the physical domain corresponding to 


te= s(t) Pest (4.52) 


with «,(7) being given by (4.51). 


time steps in the logical domain, can be recovered from (4.38) by integrating 


with respect to T, with the q’ physical time point being given by 


dr 
=¢; a: ee = Ga (4.53) 


where #(0) = 0 (time starts when particles are first emitted from the plasma 


boundary). 


Before numerically approximating the algorithm listed above by discretising 
the logical solution domain as described (using finite differences to represent 
derivatives), it is useful to list some standard derivative difference approxi- 


mations. 


e Second order difference approximations for the derivative 0 ff /0E, 


at the nodal point (jAé,kAr) for some function f(€, 7): 


— Second order forward difference approximation 


afr ee —3fe+4 a 


1 x ye Aé? 4.54 
¥ sp + O(e (4.54a) 
— Second order central difference approximation 
Of f he = et 
SS A&? 4.54b 
pe» + (ae?) (4.54b) 


— Second order backward difference approximation 


Off  3ff—4fh.it fhe 
BE i =e + O(AE?) (4.54c) 


e Second order central difference approximation for the derivative 


O fF /OE?, at the nodal point (jAé,kAr) for some function f(€, 7): 


Oi Ste ey ee 2 
3e i 2 Es + O(AE?). (4.55) 


4.3.2.3. Detailed Numerical Algorithm 


We now list a more detailed description of the steps in the implementation 


of the numerical time-stepping algorithm described above. 
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Step 1 


e Set the initial physical region size by defining the initial plasma bound- 
ary location s! = s(0) = so (noting that §1 = 0) and the target location 


LT. 


e Set the scaling parameters ¢,, Pm, and L; consequently define the 


scaling parameters v,, Eo and typ. 


e Create the initial solution vectors {f}} = {6;} = {0}, where j € (1,n+ 


1] and j € Z. 


Step 2 


e Set {of} = 0;, where k € [1,m+1] and k € Z (with m being the desired 


number of time steps). 


e Set {6%} = p,(t*), where k € [1,m+1] and k € Z, and where 


pe(t) = (1 + tanh(b(tmt — c))) 


h 
= Gp wees’ (4.56) 


for a ramped (in time) charge density at the plasma boundary, from 
(4.36f). In (4.56), the constants h, b, and c are parameters determining 


the shape of the ramp (see figure 4.3); in this case, h = 1 so that 
ps(t) € (0, 1). 
We choose a ramped (in time) charge density profile to represent the 


real physical situation, where the charge density is switched on, having 


a finite switch-on time (or ramp). 
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Figure 4.3: Function used to represent time-dependent charge density ramp. The 
parameters h, b, and c in (4.56) are indicated. The subscript n in the axis titles 


indicates a non-dimensional variable. 

Step 3 
e Choose the form of f(€,7) in (4.37). 
e Choose the form of g(7) in (4.38). 


e Select the parameters jz and v in (4.37) and (4.38). 


Step 4 
e Determine { ff } as follows: 
— When j = 1, calculate ft from (4.544); 
— When1l1 <j <(n+1), calculate fé from (4.54b); 
— When j = (n+1), calculate fé from (4.54c). 
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e Determine g* as follows: 


—Ifk=1, set g* =0; 


— Ifk > 1, calculate g* using the first order difference approximation 
i 

Step 5 
e Using {fé} from Step 4, calculate «* (Kk, at the current time step) 


from (4.51) making use of a simple trapezium rule approximation. For 


example, for some function O(é ,7), then 


fj+1 | AE fs s 
[ence SG 8), 
where pf = O((j-DAE, (k—1)Ar). 


e From (4.37), and using the value «* just calculated with {fé} from 
Step 4, determine the mapping function vector {MW oy for the current 


ime step. 
(k*) time step 


e From (4.38), and using g* from Step 4, determine the scalar time step 


controlling function value G* for the current (k“") time step. 


Step 6 


e Recover the dimensionless physical nodal locations {%;} from (4.52), 
using { fr } and «* from Steps 4 and 5. The simple trapezium rule 
approximation shown in Step 5 is used for evaluating the integral in 
(4.52). The dimensional, physical nodal locations are determined from 


{¢;} by multiplying by the constant L, established in Step 1. 
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e Calculate the vector of physical nodal speeds {ar}. = 0x/Or from the 


current and previous physical nodal locations. 


— Ifk =1 then {x*} = {0}. 
— Ifk > 1 then calculate nodal speeds using the first order difference 


approximation 


e From (4.53), recover the current, dimensionless physical time ¢(r*), 
using the vector of previous time step controlling function values {G%}, 
where q € [1,k], g € Z. The dimensional, physical time is determined 


from t(r”) by multiplying by the value t,, calculated in Step 1. 


Step 7 Using {MF} from Step 5, solve 


1 Pd" 1 AM Ad" 
Te OEGe he MEO ice ce (4.57) 
M? 0€2 = =M3 0€ O€ 
in logical space, where p* is the current charge density solution, and where 


a, is given by (4.32). To solve (4.57), we use the boundary conditions 


for Vk € [1,m+1], where ¢p and ¢g are the given conditions (4.36a) and 
(4.36b), noting (4.49). Since ¢ is known at the boundaries of the logical 
domain, we need only determine the solution to (4.57) at interior nodal points 
(i.e. for 1 < 7 < (n+1)). Thus we need only use the second order central 
difference approximations (4.54b) and (4.55) for the derivatives in (4.57). 
Using these approximations we arrive at the following system of equations 


for ok 
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e When 7 = 2, 


2 


A L Ms b 
k\2 = 5 +) TymeAg2 GNA 
(MpPAR”* * | (APARNA 
~ ot {graye + tty | 


e When 2<j <n, 
ee ee ee ee 
(Mf)?Ag? — 2AE(MP)S f°" (MP)PAg?” 
1 Mj, tk k 
$4) es — ss Oh = Oy. 
[arts 2AE(M#)3 [°F : 
e When j =n, 
Ne [| , 
VRE QREMEE Se Eee CEN ee 
1 ME Z 
ee ee Te 
fe sie} 


Step 8 Using {7é }, {MF}, and «f from Steps 4 and 5, calculate the electric 
field {Er} in logical space from (4.1d) and (4.39), by determining 


Bre =a (4.58) 
As in Step 4, the three cases 
e When j =1 calculate the derivative aG* /0E using (4.544). 
e When1l<j<n+1 calculate the derivative OG* /OE using (4.54b). 
e When j =n +1 calculate the derivative OG* /OE using (4.54c). 


apply for calculating the derivative in (4.58). 
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Step 9 Calculate the updated velocity field fort} at the next time step 
from (4.46), using {é*}, with {x}. {M}}, and G* from Steps 5 and 6. 
By rearranging (4.45), an expression for the updated velocity at the node 
(jAé,kAr) can be found using a forward difference approximation for the 


time derivative; this gives rise to the following: 


e Using a first order forward difference approximation for the time deriva- 


tive in (4.46) we have 


ER 41. @k Oo" 
ak+1 ak k Aki 
Of" = Of + ArG {+ a(e-#)% a 
J 


for the updated velocity field {oe where a is the speed of the j'® 


physical node at the k‘" time step, where 


Nk Nk Nk 


J ~~ Jj 


dE Aé 


for calculating the derivative in (4.59) also applies. 


Step 10 Calculate the updated charge density field {et at the next time 
step from (4.48), using the most recent velocity field, {ort} from Step 9, with 
{M#}, {xk} and G* from Steps 5 and 6. Rearranging (4.48), an expression 
for the updated charge density at the node (jA€, (k + 1)Ar) can be found, 
again using a forward difference approximation for the time derivative; this 


gives rise to the following: 


e Using a first order forward difference approximation for the time deriva- 


tive in (4.48) we have 


Sigs ae AOA Faas EOP Ok 
pst! = pk — H (Ga = aaa + pF 3 (4.60) 
j 
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for the updated charge density field oamee where 


ak ak ak 
Of; a) PG PFR1 


a Aé 


for calculating the derivative in (4.60) also applies. 


Step 11 


e Update the plasma boundary location by examining the electric field 
at the boundary E(s*,t*) = E* from Step 8, where s* = 8(t*); the 


electric field gradient at the boundary 
= p(t"), (4.61) 


from (4.1c) and Step 2, and the electric field curvature at the boundary 
OE;(s*,t*) 1 Ogk 
Oz ME O€’ 
from (4.61) and (4.39), where the derivative Opk /OE is calculated using 
the forward difference (4.54a). 


An expression for the size of the boundary displacement at the time 
t* can be found by determining a perturbation, As(t*) = As*, to the 


boundary using the Taylor expansion 


E(s* + As*, t*) = E(s*, t*) + As E;(3*, *) 


Ask)? ~ S 
+ OS) fisa( i) + O(A8). 


Truncating the expansion at the second order term, and making the as- 


sumption that the perturbation AS* moves the boundary to its correct 
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location (i.e. where E(8(t") + As*,#*) = 0 from the condition (4.2c)), 


gives the expression 


— Bz (5*, t*) + y/ Ez(8*, t*)? — 2E(5*, t®) Bez (5 
E33(8*, t*) 

for As*. If the system is allowed to settle to an equilibrium state of 

particle flow*, we would expect that as k — oo the plasma boundary 

will converge to the constant location §(t°) = 8, and hence also expect 

that As* — 0 as k > oo. Additionally, from the condition (4.2c), the 

term 2E(8*, t*) Fzz(5",t*) — 0 as k > oo in (4.62), indicating that we 


must take the positive square root (or As* + 0 as k — oo). 


— If (E;(8*, t*)? — 2F(5*, t*) Esa (5*, t*)) < 0 then set AS* = 0. 


Step 12 


e Recover the dimensional solution variables from the dimensionless def- 


initions (4.27a) to (4.27f). 


e Increment & and return to Step 4. 


The Boundary Perturbation Stabilisation Parameter. After imple- 
menting the algorithm, it was observed that execution would proceed to 
relatively long times only if the calculated boundary perturbation at each 
time-step was small. Therefore, to prevent the algorithm failing and thus 
obtain analysable results, an additional parameter, 6, was introduced, which 
pre-multiplies the calculated boundary perturbation, thus reducing its mag- 
nitude at each step; this allowed the numerical algorithm to proceed unhin- 


“As in the time independent case. 
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dered, but with the unwanted payoff of the predicted boundary movement 
being un-physical. 

To update the boundary location at each time-step, the calculated bound- 
ary perturbation, As*, (given by (4.62)) is multiplied by the small parameter 


6 in that the updated boundary location at the (k+1)™ time step is given by 
gett = Bh + 5A8*, 


The introduction of 6 in this way allows the boundary perturbation to be 
reduced to a level that allows the numerical method to proceed without fail- 
ure, although the calculated time-dependent boundary location is unlikely to 


be representative of the true physical case. 


The algorithm detailed above has been coded in MATLAB, and we now 
present results from its operation on a test case chosen, in hindsight, to 
not only represent experimental operation of a neutron tube, but to also of- 
fer physically representative results within a reasonable computational time 
frame. In this test case, the charge density at the plasma boundary is ramped 
to a maximum over a period of tens of nanoseconds. However, prior to de- 
tailing these results, we derive a CFL condition ([26], [24]) applicable to the 
scheme employed; this ensures stable choices of logical spatial and time step 


size. 


4.3.2.4 CFL condition 


Implementation of the algorithm developed in Section §4.3.1 required exper- 
imentation, as one would expect. The logical spatial and temporal step size 


must be chosen to avoid numerical instability caused by violation of the CFL 
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condition ({26], p87) appropriate to the first order, upwinded difference ap- 
proximations, (4.59) and (4.60). To determine the condition, and due to the 
dynamic nature of nodal point locations in the moving physical region, an 
explicit expression for it is now evaluated on the logical, fixed region, being 


initially given by 
Ag 


< — : 
Se (4.63) 


Ow 
where AE is the logical nodal spacing and Az is the logical time step. The 
wave-speed %,, is deemed to be the maximum speed of the initial particle 
wavefront as it crosses the logical region, and can be determined by consid- 
ering the maximum wavefront speed in the non-dimensional physical region, 
as follows. 

Particles initially released from the plasma boundary form a wavefront, 
which undergoes uniform acceleration across the non-dimensional physical 
region, caused by the initial charge free uniform electric field . Using the 
familiar expression for the distance travelled by a particle under uniform 
acceleration, the time ¢, taken for the wavefront to traverse this region is 


found by determining the positive root of 
1 = O;t. + =at? (4.64) 


where 0; is the initial non-dimensional particle speed, and a is the non- 
dimensional uniform acceleration of the wavefront, with the 1 on the left 


hand side of (4.64) being the region size. Thus 


ee 
t= (-# + 4/0? + 21 (4.65) 
a 


The maximum wavefront speed in the non-dimensional physical region is 


1Ae 


therefore given by 


Oy = Gt. +0; 
= 4/0? + 24, (4.66) 


from (4.65), where the non-dimensional acceleration @ can be determined 


from (4.35b), since 


G2 sO OU 
di 0% «Ot 


_E 
me 


By applying the boundary conditions (4.36b) and (4.36c), the initial uniform 


electric field strength formed within the charge free non-dimensional region 
is 


Eo = ($s — br)/(&r — 5(0)) 
=e 


so that the maximum wave-speed (4.66) can be written 


Ow = 4/0? +1. 


(4.67) 


The mapping of physical non-dimensional time ¢ to logical time 7, using the 


controlling function (4.38), indicates that the maximum rate of change of t 
with 7 is given by 


OT 
or by 


(4.68) 
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If the time taken for the wavefront to cross the non-dimensional region is 
given by 

it 
Las Se. (4.69) 


then the minimum time taken for the wavefront to cross the logical region 


will be 
“ 1 
1 = 
Vw 
Tes 
= eee 4.70 
7 (4.70) 


from (4.68), so that the maximum wave-speed in the logical region is 


dy =n/2 +1 (4.71) 


from (4.67), (4.69) and (4.70). Additionally, 


VU; FS 


lI 
ie 
8/8 
3 

S 


from (4.28) and (4.36d), and when combined with (4.71), this gives rise to 


= Be Os 
“Y 2abm Ar’ 


from (4.63); this is applied to the logical spatial and time step sizes. 


the condition 


4.3.2.5 Test Case 


To test the algorithm, the initial conditions listed at the start of §4.3, above, 
were applied. Here, the potential difference applied across the solution region 
was held constant in time, with particles emerging with a constant initial 


velocity (or energy); this was to represent an experimental setup, although 
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numerically, these conditions could have been time-dependent. The initial 
region size was adjusted in successive runs to accelerate convergence (towards 
the steady state solution) as numerical experiments progressed. Additionally, 
the maximum charge density (pm = Jo/v;) at the plasma boundary was 


selected to represent a typical value within an actual device, and the functions 


F(E,7) = pl§,7) (4.72a) 
g(r) = a(é(r)), (4.72b) 


were chosen as the governing functions in the mapping function, (4.37), and 


time step controlling function, (4.38), respectively. 


Source Charge Density Ramp. In this case, as a representation of an 
experimental increase in charge density at the plasma boundary, the applied 


ramp is increased from zero to the constant value 


ps(t) = Jo/vi 
— Pm 
= const (4.73) 


with the increase occurring over a period of time of the order of tens of 
nanoseconds”. Within the algorithm, this is effected using the function (4.56) 


shown in Figure 4.3, where the non-dimensional charge density is increased 

5In comparison with experiment, this is actually at the short end of the ramp time 
frame, but was chosen this way after the necessary inclusion of the parameter, 6. The 
inclusion of 6 causes the numerical method to execute very slowly and so to be able to 
examine the effects of a ramped increase in charge density reminiscent of any kind of 


experiment, the shortest applicable ramp time frame was chosen. 
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from zero to one over a period of non-dimensional time of maximum order 
ten. In this case (by referring to (4.56) and Figure 4.3) the constants }, c, 


and h are chosen to be 


b=1x 10's 1, 
3 
C= > 
b 
=3x108s 
h = Pm, 


where pm is given by (4.73). These parameters cause a ramp rise time of 
~ 6.88 non-dimensional units, or 20 ns; the value t,,, = 2.91 ns follows from 


(4.28) and (4.29) (with L = 0.009m from table 4.1). 


4.3.2.6 Results 


Application of the physical conditions listed in §4.3 above, including the 
short-timescale charge density ramp (detailed above) along with a logical 
time-step size of 5.5 x 1073 non-dimensional units (corresponding to ~ 182 
non-dimensional time steps per unit ¢,,, or an initial physical time step of 
1.6 x 1071's), yielded results with a number of differing features; some of 
these are explained in the following paragraphs. Results are displayed in 
physical space. 

It is pointed out here, that in order to obtain a suitable set of results, a 
number of numerical experiments were performed and these resulted in the 
set of parameters listed in Table 4.1 that allowed the algorithm to execute for 
> 20000 time steps. Here, the time step size shown was chosen for stability, 


and allowed the procedure to proceed without failure to a final physical time 
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Number of time steps. 


60 Number of spatial steps. 


5.42 x 1073 units | Logical time step size. 


0.015 Parameter to adjust adaptive phys- 


ical nodal density (see (4.37)). 


1 x 107° Parameter to adjust adaptive phys- 


ical time step size (see (4.38)). 


Parameter to adjust physical time 
step size outside of adaptive regions 


(see (4.38)). 


Parameter to adjust amount of 
boundary movement at each time 


step. 


Region start size - final run. 


Table 4.1: Numerical parameters for the one-dimensional time-dependent numer- 


ical method 1 test problem. 


of 320 ns; at this time the applied source charge density ramp has stabilised. 
Numerical experiments showed that the number of nodes chosen (n = 60) 
was a compromise between speed and accuracy, with higher numbers of nodes 
significantly increasing execution time and significantly fewer nodes (<30) 
causing a violation of the CFL condition, (4.72), or other instability; these 


results are not detailed here. The initial region size of L = 0.009 m is typical 
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of a neutron tube source shield to target separation distance. The parameter 
c, was set to a default value of 1. The parameters ys and vy in (4.37) and 
(4.38) were adjusted to aid speed of execution with the particular set of 
initial conditions used over a number of experiments. The resulting values 
chosen caused the most efficient execution of the algorithm in conjunction 
with the other parameters listed in Table 4.1. 

The following Figures detail features from the calculated solution, and in- 
clude comparisons between the time-dependent solution and known analytic 


solutions. 


Figure 4.4 shows the calculated surface p(x,t). The ramped increase in 
charge density at the plasma boundary is visible, as is the adaptive time 
stepping procedure, with widely spaced time steps to the left of the surface 
being due to zero boundary movement over this region of time; this is caused 
by the discriminant in (4.62) being negative, invoking the error trap condi- 
tion detailed in Step 11 of the difference algorithm, above. Figure 4.5 shows 
level curves of charge density over the first ~ 93 ns. The plasma boundary 
is identified by the solid black line to the left of the Figure, and the time at 
which boundary movement begins (~ 33ns) is clearly visible®. At this point, 
the sudden, rapid boundary movement is propagated within p(z,t), and this 
is likely to be due to sudden changes in the time step controlling function G 
(which is dependent upon the rate of change of boundary position, (4.72b) 
and (4.38)). A reduction in the logical time step size would reduce this effect. 
The separation characteristic separating the charge free region (ahead of the 


°It is pointed out that this time is largely a numerical artefact of the method employed 


and not necessarily a reflection of the true, physical situation. 
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Figure 4.4: The charge density surface p(x,t). Output is shown at every 200 time 


steps. 


advancing wavefront) from the charge infused region behind it, is superim- 
posed. The entire solution region is charge infused within ~ 6 ns, although 
due to the ramped initial condition at the plasma boundary, the charge den- 
sity within the region is relatively low at this time. To identify a level curve 
close to the separation characteristic, a total of twenty logarithmically spaced 


levels are displayed in the range 10~*° < p(z,t) < 10°, where 


b = log max [p(z,t)] } , 


x,0<t<tmax 
and tmax = 93ns. This is further demonstrated by examining the solution 


up to the earlier time of ~ 10 ns in more detail, shown in Figure 4.6. Again, 
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similarly calculated level curves of charge density are displayed, where the 
curve towards the base of the Figure follows the superimposed separation 


characteristic closely. Figure 4.7, shows solutions extracted from the electric 


x 10° p(x,t) Level Curves 
T T T £ T T 
9P | Charge Density 1 
om Separating Characteristic 
gb Plasma Boundary I 


Time (s) 


Distance (m) <8 


Figure 4.5: Logarithmically spaced level curves of constant charge density over 


the first ~ 93 ns. 


scalar potential (voltage) surface ¢(x,t) at different times. The progression 
from the initial charge free voltage solution (at 0 ns) towards the curved 
steady state solution is apparent, with the 93.08 ns curve appearing to be 
very close to the analytic (black “dash-dot”) steady state solution. Finally, 
Figure 4.8 shows the variation in plasma boundary location s(t) with time. 


The red curve is the calculated plasma boundary location, whereas the green 
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Figure 4.6: Logarithmically spaced level curves of constant charge density over 


the first ~ 10 ns. 


“dash-dot” curve is the calculated plasma boundary perturbation As(t); this 
function has been multiplied by the reciprocal of the parameter 6 to increase 
its magnitude to enable display within the same graph as s(t). By setting the 
initial region size to 0.09m, the numerical algorithm has taken ~ 93ns for the 
calculated boundary position to reach the analytic steady state location; the 
algorithm was terminated at this point. Clearly, from the gradient of s(t) as it 
reaches the steady-state boundary location, it is apparent that the calculated 
plasma boundary does not settle at this location and will overshoot somewhat 


as time progresses. By observing that after about ~ 50 ns, the perturbation 
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Figure 4.7: The variation of electric scalar potential (voltage) with distance across 
the region at different solutions times. The analytic steady state solution and 


steady boundary location are shown. 


As(t) is approximately linear, an estimate of when the boundary will begin 
to recede again can be made. Performing a linear regression on the linear 
portion of the function As(t) allows the point at which perturbations become 
negative to be estimated; this is calculated to be ~ 219.6 ns, as shown in 
Figure 4.9. However, the increased electric field magnitude throughout the 
region, caused by the boundary movement and consequent reduced region 
size, may cause the boundary to begin to settle at an earlier time than this. 


Having observed calculation times for differing numbers of time-steps (Figure 
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Figure 4.8: Variation in plasma boundary location s(t), and plasma 


perturbation As(t), with time. 


boundary 


4.10), it is estimated that to reach a physical solution time of ~ 219.6 ns 


using the algorithm (in its current form) on the fastest desktop computer 


available, execution times will be ~ 8 hours, taking ~ 35000 time-steps 


(Figure 4.11). Due to arrangements surrounding the use of the computer, 


execution times of this size have not been possible. However, reaching this 


point in the calculated solution is largely academic as a consequence of the 


introduction of 6; the predicted boundary movement is unlikely 


accurate representation of the actual physical process. 
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Figure 4.9: Estimated time at which the boundary perturbation As(t) becomes 


negative. 


Adaptive Nodal Density Procedure. ‘To test the adaptive nodal den- 
sity procedure, a number of runs were conducted for each case with differing 
values of the parameter py. If js was set to be of order ten, or greater, prob- 
lems in calculating derivatives at the plasma boundary were observed. An 
example of this is shown in Figure 4.12, where the electric field (in physical 
space), derived from the calculated potential field via (4.58), has an apparent 
discontinuity at the plasma boundary. The failure is a consequence of the 
mapping function WM ; (from (4.37)) becoming relatively small at the bound- 
ary and is caused by the gradient of the controlling function (4.72a) becoming 
very large there. 


However, the relatively gentle increase in charge density (and hence con- 
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Figure 4.10: Estimated algorithm execution time to reach the point at which the 
boundary perturbation As(t) becomes negative. Points marked “observation” were 


taken from the current Case. 


trolling function) at the boundary allows the procedure to work without 
failure, and setting 4s = 10 causes the nodal density near the plasma bound- 
ary, and hence the effect of the calculated derivative pe, to be high. This can 
be seen in Figure 4.13, where the calculated physical potential, ¢(x,t), across 
the acceleration region is shown at three different solution times. Nodal loca- 
tions are highlighted by markers situated along each curve, and at the later 


time, these points coalesce towards the plasma boundary location to the left 
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Figure 4.11: Estimated number of time-steps to reach the point at which the 
boundary perturbation As(t) becomes negative. Points marked “observation” were 


taken from the current Case. 


of the graph; this is a consequence of the spatial adaptive nodal density 
procedure. 

A further demonstration of the adaptive procedure can be seen in Figures 
4.14 and 4.15. Here, the locations of the first fifteen nodes, in physical 
space, are displayed as a function of time and the variation in physical nodal 


separation is displayed at varying solution times. In both graphs, nodes are 
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Figure 4.12: The electric field surface E(x,t), with adaptive parameter = 100 


(region start size of L = 0.009 m). Output is shown at every 40 time steps. 


initially equally spaced, seen to the left of Figure 4.14 and by the blue line in 
Figure 4.15, but as time progresses they become rapidly concentrated towards 
the plasma boundary location, again highlighting the adaptive nodal density 
controlling. In Figure 4.15, this is manifested at later times by the nodal 
separation rapidly decreasing towards node 1 at the left. Additionally, away 
from the plasma boundary at the left, the nodal spacing remains constant 
as a consequence of the function fz(€,7), being effectively zero there (this 
is f(€,7) in (4.37)). The green (3.6ns) curve in Figure 4.15 shows some 
instability, although this decreases with time and has disappeared after a 


solution time of about 6ns. 


132 


x 10° Variation of Potential with Distance with Nodal Locations (1=10) 
12 T T T T T T T 


—e0ns 
—#- 5.97 ns 
—7— 46.62 ns 


Potential (V) 


0 if Ll i Ll L L Ll i 


4 5 
Physical Distance (m) 


Figure 4.13: Variation in potential ¢(z,t), with x at three specific times. Setting 
the adaptive parameter, jz, to be equal to ten leads to increased nodal density 


towards the left of the curves, as time increases. 


4.3.2.7 Method 1 - Summary and Conclusion 


Whilst the algorithm presented here offers an apparently viable solution, in 
that the advancing wavefront follows the analytic separation characteristic 
closely (See Figure 4.6, for example), in order to offer stability to the method 
the small, multiplicative parameter, 6, must be introduced into the plasma 
boundary calculation thus rendering the predicted boundary movement un- 
physical. Without such a parameter, the boundary movement calculated by 
(4.62) causes rapid failure of the numerical method, with this being due to 
the presence of £;;(8",t") on the denominator of (4.62). Particularly at 


early solution times, the calculated potential, d(@), across the solution re- 
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Figure 4.14: Variation in nodal location for the first fifteen nodes in physical space. 
Nodes rapidly move from an equally spaced distribution to being concentrated near 


the plasma boundary (where the charge density gradient is greatest). 


gion has little curvature. Consequently, the term Eee’. t*), corresponding 
to the third derivative of the solution potential at the plasma boundary, is 
very small, in turn causing As’ to be large. Such large changes in boundary 
location cause the solution region to become distended, with the consequent 
effect being that the mapping function, M, becomes large rendering the di- 
mensionless, logical version of Gauss’ law, (4.44), imbalanced at successive 
time steps; this is manifested as a rapid onset of instability within the nu- 
merical procedure, causing it to fail. 

A further consequence of the necessity of 6 is that the algorithm executes 


very slowly. Indeed, in order to examine the effects of an applied source 
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Figure 4.15: The variation in nodal separation with node number at different 
solution times. The ordinate value at node 2, for instance, refers to the distance 


ck — * (for time-step k). 


charge density ramp that can be replicated at all by experiment, a relatively 
rapidly increasing ramp must be applied. Whilst this is (just about) experi- 
mentally viable, it is not typical of a normal neutron tube operation. 

In order to avoid the issues within the first method, a new approach is now 
proposed, where an expression for the rate of change of boundary movement 
is determined in an entirely different way. This new approach does not suffer 
the same problems experienced here, and thus allows a more experimentally 


representative charge density profile to be applied within the model. 
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4.3.3. Numerical Method 2 for the Determination of the 


Time-dependent Solution 


As mentioned above, a new approach to the time dependent problem has been 
developed to avoid the problems experienced in the original method; this is 
now described along with an associated time-stepping algorithm. We begin 
by first manipulating the original one-dimensional, dimensionless system of 
equations, (4.35), to reduce their number by one and then further manipulate 
the system to yield an analytic integral expression for the rate of change of 
dimensionless plasma boundary location, §’(t). The approach is described as 


follows, beginning by first recalling (4.30) and (4.35a), which are 


2 dd 
E=-— 
Ox’ 
and 
fag 
dg? 
respectively. 


By noting (4.30), we can integrate (4.35a) to observe that 


BGA = orf pg,f) a, (4.74) 
3(8) 
using the condition, (4.36c). Then, substitution of (4.74) into (4.35b) gives 


rise to the coupled pair of equations, 


Oo 20g fg 0 

Res Shee dg — 05m (4.75a) 
dp __ api) 

oe (4.75b) 


for the dependent variables f(Z,t) and i(Z,t), with the boundary condition, 
(4.36c), being used in (4.74). 
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An expression for the rate of change of plasma boundary location, $’(t), 
can also be determined by initially integrating (4.74), over the solution do- 


main, giving 


[,* dé = or — os 


= -a, if : (fr — @)p(&, t) dé, (4.76) 


by the replacement lemma ([38], p6). 


Taking the time derivative of (4.76), we find 


by Leibnitz Integral Rule ([36]), and since 
a(pe) _ a 
0 Ot’ 
from (4.35c), the integral on the right hand side of (4.77) can be expanded, 


giving 


= (&r — 5(t))p53"(t). (4.78) 


Evaluating the boundary term on the right of (4.78), we arrive at the equation 
(i) 1 "6% (4.79) 
ps(@r — 8(t)) Jag 


for the rate of change of dimensionless plasma boundary location at f. 
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4.3.3.1 Detailed Numerical Algorithm 


The equations, (4.75), in addition to (4.79) can now be used in a time step- 
ping algorithm for the determination of the solution of the time-dependent 
plasma boundary problem in one dimension. The new proposed algorithm is 


as follows: 


Step 1 The initial solution region is set with the initial plasma boundary 


location being defined as §(0) = So; 


Step 2 Set 6(%,0) = 0(%,0) = 0, where & € (5, 7]; 


Step 3 Set p(s,t) = ,(t) and 0(S,t) = %;, from (4.36d) and (4.36f), where 
fs(t) is the dimensionless time-dependent function (in this case given by 
the ramped function, (4.56)) representing the time-dependent variation in 
current density at the plasma boundary during tube operation, and v; the 
initial dimensionless particle velocity (a constant in this case, although not 


necessarily so); 


Step 4 Calculate the updated velocity field, 6(Z,t**1), from (4.75a), us- 
ing the current velocity field, 6(Z, t*), and the current charge density field, 


p(z,t*). This is done by discretising (4.75a) as 


ght! — pha Ad otf” jae a ah 
ey DA i oe | f° 


where the integral is approximated by the trapezium rule, 


a} Agk = 
[aes a area) rsa 


1=2 
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with Az* being the nodal spacing at t*. Here, the spatial derivative aur /O£ 


is approximated by the first order difference 


Step 5 Calculate the updated charge density field, p(, t*+'), from (4.75b), 
using f(Z,t”), and 0(z, t**') from Step 4. This is done by discretising (4.75b) 


as 
ke k+1 
j OX 
Here, the spatial derivative 0 (2 pro - /O£ is approximated by the first order 
difference 
xk k ~k wk 
OBFOF Pe — Faia 
Ox Dec 


Step 6 Determine the rate of change of boundary movement at the current 


time-step, t*. This is done by evaluating 


dgk Au Ze (pro; + pk ORT) a pa 2Pj apn} 
i.) - —— - Gee -— > 


: (4.80) 
from (4.79), where 5* = &(t*). Here, the integral on the right hand side of 
(4.79) has been approximated in (4.80) by the trapezium rule, as in Step 4. 
An approximation to the boundary location at the (k+1)™ time step can 
then be written 

dst 


xk+1 ~ Aim. 
S gy de’ 


Step 7 Return to step (4) and repeat until the solutions in successive time 


steps are identical to within a small tolerance. 
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4.3.3.2 Test Case 


To test the new method, the physical conditions listed at the start of §4.3, 
above, were applied. As in the previous method the potential difference 
applied across the solution region was held constant in time, with particles 
emerging with a constant initial velocity (or energy). As a consequence of 
a significant increase in computational speed, a source charge density ramp 
that is more typical of experiment could be used and this is detailed as 


follows. 


Source Charge Density Ramp. Owing to long algorithm execution times, 
in order to have any realistic representation of experiment the previous model 
used a source charge density ramp whose rate of increase was rapid enough 
that the effects of the ramp could be captured within a reasonable compu- 
tational time frame. Such a rapidly increasing source charge density is at 
the short end of the experimental time frame, and is not typical of neutron 
tube operation. The new method presented here overcomes difficulties expe- 
rienced in the previous case and so to test it, a more typically representative 
charge density ramp is applied here. Again, referring to the function (4.56) 


shown in Figure 4.3, the constants 6, c, and h are chosen to be 
b= 5Al< 10°%s, 
6=2.%10-%s 
h = Pm, 


where p,, is given by (4.73). These parameters cause a ramp rise time of 
~ 114.4 non-dimensional units, or 370 ns; the value t,, = 3.232 ns follows 


from (4.28) and (4.29) (with L = 0.01m from Table 4.2). 
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4.3.3.3 Results 


The application of the physical conditions, including the more representative, 
longer timescale charge density ramp in addition to a logical time-step size of 
4.64 x 10~? non-dimensional units (corresponding to ~ 215 non-dimensional 
time steps per unit t,,, or an initial physical time step of 1.5 x 107! s) and 
a total time of 900ns, yielded results with a number of differing features; 
some of the more salient features are explained in the following paragraphs. 
Other parameters used are listed in Table 4.2, below, where it is noted that 
the adaptive features incorporated in method 1 are not implemented here. 
The algorithm executed to completion within a time frame that was of the 


order of minutes, rather than hours (as in the original algorithm). The pa- 


Number of time steps. 


75 Initial number of spatial steps. 


4.64 x 107? units | Logical time step size. 


0.01 m Region start size - final run. 


900 x 107° s Calculation termination physical 


time. 


Table 4.2: Numerical parameters for the one-dimensional time-dependent numer- 


ical method 2 test problem. 


rameters chosen allow a solution showing convergence towards the analytic 
steady-state solution to be generated, with the calculated location of the 


plasma boundary clearly asymptoting towards the steady state location; the 
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accuracy of the final solution increases with increasing numbers of spatial 
nodes. The time step size was chosen to facilitate this for differing numbers 
of spatial nodes without violating the CFL condition, (4.63) (which is also 
applicable in this case). Initially, 75 spatial nodes were used, with this num- 
ber increasing in steps of 25 to 150 nodes. At the final solution time of 0.9ys, 
a typical neutron tube source will have undergone switch-on with the source 


charge density being approximately constant at this time. 


The following Figures detail features from the calculated solution, where 


results are displayed in physical space. 


Figure 4.16 shows the surface p(x,t), where the charge density ramp, in- 
creasing from right to left, is clearly visible. There is no obvious instability. 
Figure 4.17 shows level curves of charge density for the entire calculated 
solution. The plasma boundary is identified by the solid black line to the 
left of the Figure, where it can be seen that boundary movement begins 
immediately. The entire solution region is charge infused within ~ 7 ns. 
Figure 4.18 shows the surface (x,t) (extracted from p(%,t) via (4.35a) 
and re-dimensionalised). The change in plasma boundary location (at the 
top of the surface) with time (increasing from left to right) is visible, as is 
the gradual change in curvature of the potential across the solution region 
from the initial, linear solution at t =0 (to the left of the surface) to the 
final curved solution at t = 900ns (at the right of the surface). To examine 
the final solution, it is displayed in Figure 4.19 as the solid black line; it is 


surrounded by red circular markers, which represent the analytic steady-state 
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Figure 4.16: The charge density surface p(x,t). Output is shown at every 600 


time steps. 


solution (extracted from (3.28)). The two curves (black line and markers) 
appear very close together, and so to further investigate the accuracy of 
the final time-dependent solution, a number of calculations were performed 
with differing numbers of nodes. For each calculation, the final calculated 
time-dependent potential solution (at t = 900ns) was compared to the exact 
steady-state solution calculated at the same nodal point locations (again 
using (3.28)). The error in calculated final solution (as a percentage of the 
maximum potential value, ¢,,) for each of these calculations is shown in 
Figure 4.20. Here, it can be seen that the calculated solution error decreases 
with increasing numbers of nodes. It can also be seen that the largest error 


in solution occurs towards the plasma boundary (to the left of the graph) 
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Figure 4.17: Equally spaced level curves of constant charge density for the entire 


solution (up to 0.9ps). 


where the solution curvature is at its greatest; this is to be expected and 
would be reduced if an adaptive scheme, such as that used in method 1, were 
used here. As expected, the error in the calculated solution is zero at the 
fixed target location to the right of Figure 4.20. 

Finally, Figures 4.21 and 4.22 show the calculated boundary location and 
perturbation as functions of time. Figure 4.21 clearly shows the time depen- 
dent boundary, s(t), asymptoting to the analytic steady-state plasma bound- 
ary location, with Figure 4.22 showing the boundary perturbation asymptot- 


ing to zero with time. 
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Figure 4.18: The potential surface p(x,t). Output is shown at every 600 time 
steps. 
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Figure 4.19: The analytic steady-state solution (red circles) with the calculated 


potential solution (x,t) at t = 900ns (black solid line with 150 nodes). 
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Error in Calculated Solution Potential with Distance at 900ns 
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Figure 4.20: Error in final calculated solution potential as a function of distance. 


dm 


The error is calculated as error% = , where ft, is the final time 


(900ns in this case) and ¢,(x) is the analytic steady-state solution. 


4.3.3.4 Method 2 - Summary and Conclusion 


The new approach detailed here, for the solution of the one-dimensional time- 
dependent problem, avoids solving Gauss’ Law, (4.1c), numerically, which, 
in combination with the method used for the determination of the plasma 
boundary perturbation, is believed to be the source of the instability ob- 
served in the original method. Furthermore, a more robust, analytically 
based method for determining the plasma boundary location has been devel- 
oped. 

The new approach has been implemented as a numerical time-stepping 
scheme using MATLAB, and has been tested by applying conditions typically 


seen in a neutron tube application; this was barely possible with the original 
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Figure 4.21: The variation in boundary location, s(t), with time. The calculated 
boundary asymptotes to the analytic steady-state location as the source sharge 


density ramp reaches its maximum value. 
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Figure 4.22: The variation in boundary perturbation, As(t), with time. The per- 
turbation, calculated from (4.79), asymptotes to zero as the source charge density 


ramp reaches its maximum value. 
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method, as a consequence of slow execution times. Results clearly indicate 
that, in this case, as the source charge density ramp settles to its peak height, 
time dependent solutions converge towards the analytic steady-state solution, 
as would be expected. The calculated plasma boundary location is seen to 


clearly converge upon the steady state plasma boundary location. 


4.4 Summary and Conclusion to the One-dimensional 
Time-dependent Problem 


Several attempts have been made at solving the one-dimensional system, 
(4.1) (subject to the conditions, (4.2)), with varying degrees of success. The 
problem was initially attacked analytically, yielding a solution that is valid 
in the charge free region separated from the charge infused region by a par- 
ticular characteristic curve. Further analysis has also yielded the form of an 
analytic solution in the charge infused region, to within constants, although 
the determination of the constants appears difficult. 

A first numerical approach, whereby the system, (4.1), is transformed to 
a fixed, logical grid (as a consequence of the physical solution region being 
time-dependent), upon which it is discretised in a standard way using finite 
differences, yielded an apparently viable solution. However, in order for 
the method to proceed without instability, an unphysical scaling parameter 
was introduced, dramatically reducing the the boundary movement at each 
time-step and causing the method to proceed very slowly. It was concluded 
that the unphysical nature of the calculated boundary movement was not 


acceptable and so a new method was devised to circumvent the difficulties. 
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The new method removed the observed source of difficulty by integrat- 
ing Gauss’ Law (as in (4.74)), and incorporating it into the Lorentz force 
equation, thus reducing the system of equations by one to the pair, (4.75). 
A different method for determining the boundary movement was also deter- 
mined, and this was done by again integrating Gauss’ Law. Manipulation of 
this, whilst utilising the continuity equation, (4.75b), yields the expression, 
(4.79), for the boundary velocity; the boundary movement at each time-step 


is determined from this expression. 


We now move on to the more difficult two-dimensional case, where a method 
is proposed for the determination of the two-dimensional steady state solu- 
tion. This method is based upon the more successful second one-dimensional 
numerical method, and splits the Lorentz, vector equation, (2.17a), into its 


two components, solving each one separately in alternating directions. 
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Chapter 5 


A Planar ‘Two-dimensional 


Steady-state Solution 


One of the key assumptions within the study of the one-dimensional problem 
was that, within the acceleration region of a neutron tube, particle acceler- 
ation due to magnetic fields can be disregarded. As a result of the study 
of the one-dimensional time-dependent problem, an approach to the two- 
dimensional steady-state problem is now proposed. The proposed solution 
strategy has been implemented, giving results that are entirely plausible and 
which show convergence, with a little further work being required to perfect 
the method. 

Within this chapter, we seek a two-dimensional steady-state solution to 
the plasma boundary problem introduced in §2.3, by solving the system of 
time-dependent equations, (2.17), on the region shown in Figure 2.3 (which 
is reproduced here for convenience and labelled Figure 5.1), with time being 
used in a pseudo manner (it not representing actual time) to advance an 


iterative method. 
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It is pointed out that the strategy adopted in this chapter calculates 
the deuteron beam entirely within the solution region. This is in contrast 
to the alternative approach of solving only within the beam itself, where 
the edge of the beam (which would correspond to the edge of the domain) 
would represent a free boundary. This discarded alternative is potentially 
more difficult, because of the need to locate a second moving boundary that 


represents the edge of the beam. 
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Figure 5.1: Region of interest, R. The curved line is a representation of the edge 


of the ion beam, and separates the zero and non-zero charge density regions. 
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5.1 Two-dimensional Time-dependent System 


In order to develop the two-dimensional steady-state solution to the plasma 
boundary problem, we firstly restate the original simplified (neglecting mag- 


netic effects) system of equations as 


mio v4 oh = gE, (5.1a) 
V - (pv) = — (5.1b) 
VivVe= = (5.1¢) 
and 
E=_V¢, (5.1d) 


from (2.17) and (2.18), where t > 0, say, and where the constants g, m, and 
€9, along with the dependent and independent variables, are as defined in 
82.1 and §2.3.1. We then write this system in a two-dimensional Cartesian 
form, utilising the fact that the electric field can be written as an integral, 
from (5.1d) and (5.1c), as in §4.3.3. 

The time dependence within this system is to be used to iterate towards 
a steady-state solution on the planar region shown in Figure 5.1 (which is 
symmetric about the vertical axis). The system is solved in conjunction with 
appropriate extensions to two dimensions of the conditions (2.20) given in 


§2.3.2, namely 


o(x,t) = ds, t>0,x € Si, US» (5.2a) 
o(x,t) = or, t>0,xeES4 (5.2b) 
Vo(x,t)-n, =0, t>0,xeES; (5.2c) 
V0(x,t)- ns =0, t>0,x€Ss (5.2d) 
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p= plx t>0,xES; (5.2e) 


pHa, t > 0, x € So, (5.2f) 

p=), t=0,xER, x €S,USz, (5.2g¢) 

V p(x, t) ns = 0, t>0,x eS; (5.2h) 
v =0;(x)n, t>0,xES; (5.21) 

v=0, t=0,xER, x€¢S, (5.2j) 

v-n; =0, t>0,xeSs; (5.2k) 


where {S;} are the boundary segments enclosing the region R, shown in 
Figure 5.1. Here, S; is the (moving) plasma boundary, S2 is the top of 
the plasma cup (held at the same electric scalar potential as the plasma 
boundary, (5.2a)), and Sy, is the target (held at the constant electric scalar 
potential (5.2b)). As explained in §2.3.2, the condition (5.2c) indicates that 
the electric field normal to the plasma boundary is zero, which implies that 
the electric field at the boundary is zero, since the electric field tangential to 
the boundary surface must be zero as a consequence of (5.2a). 

The condition (5.2i) indicates that particle release speed is a function of 
spatial location only, and that particles are released in a direction normal to 
the plasma boundary surface (82.3.2), whilst the condition (5.2e) indicates 
that the charge density at the plasma boundary surface is also a function 
in spatial location only; it is zero, (5.2f), on the metal plasma cup. These 
conditions are appropriate for the steady-state case. The conditions, (5.2j) 
and (5.2g), indicate that the region is initially charge free, and the final 
conditions, (5.2d), (5.2h) and (5.2k) state that there is no change in @ and 


p across the line of symmetry, Ss, in addition to there being no velocity into 
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the boundary. 

In this case, the vector x = {x,y} is the vector of independent variables, 
where x € [0, xr], and where y € [0, yr]. The right hand boundary segment 
S3 is located at (xp, y), and the segment S, is located at (x, yr). Clearly 
XR must be sufficiently far away from S; that the beam, which is inclined 
to expand as it traverses the domain as a consequence of mutual charge 


repulsion within it, that it does not strike S3. 


5.1.1 Two-dimensional Cartesian Form 


To write the system of equations, (5.1), in a Cartesian form suitable for the 
solution strategy proposed in this chapter, we initially select the y axis to be 
parallel to the axis of symmetry in Figure 5.1, and then proceed to simplify 
the Lorentz equation, (5.la). A simplification can be achieved by setting 


p = q =v in the vector identity 
V(p-q) = (Pp: V)a+ (a: V)p + p x (V x q) +a x (V x p) 
(see [29]), and writing (after rearranging) 
(v-V)v= 5V(vev)—v x (V x v) 
= sv (v-v), (5:3) 


since the assumption has been made that there are no magnetic effects and 
hence no vorticity exists within the particle flow (i.e. V x v = 0). Substi- 


tuting (5.1d) and (5.3) into (5.1a) gives 


sv (v-v)+ (=) =-4v¢, (5.4) 


where for an iterative scheme to determine the steady-state solution, we could 
neglect the time dependence in velocity (using the time dependence in (5.1b) 


to progress the scheme), thus allowing (5.4) to be integrated to give 


= =o +4, (5.5) 


where v, and v, are the x and y components of the velocity vector, v, re- 
spectively. The constant c; can be determined by noting that ¢ = ¢s when 
|\v| = v; (on S1), from (5.2a) and (5.2i). In practice, however, (5.5) proves 
difficult to use, particularly at early stages within the calculation, since the 
extraction of v,, say, from it often gives complex results. As such, it is 


preferable to split (5.4) into its two Cartesian components, 


and 


At this point, we can determine integral expressions for the right hand sides 
of (5.6) in a similar manner to the one-dimensional formulation given in 


§4.3.3. 


By considering (5.1c) in two-dimensional Cartesian form and rearranging, 


we can write 


deettsy) =~ (FEY + bp ea)), (5.7) 
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where, by integrating (5.7) for a given yo location within the region, we can 


either write 


/ baa ee 


= P p(x, Yo) rg 
= a €0 + Oy 


where s, is the x location of the plasma boundary at yo, and where s,, is the 


) dé, 0 is Yo x Sh, 
&,YO 


maximum height of the plasma boundary above the source cup, or 


/ bea(®, Yo) AE = o2(2x, yo) 


0 €0 Oy? 


where x = 0 corresponds to the symmetry line to the left of Figure 5.1, with 


dé, Sh < Yo S UT; 
&,YO 


yr being the target height (the height of the boundary S,). Within both of 
these expressions, the value of ¢,(x, yo), evaluated at the lower integration 
limit, is zero. This is since E(s,) and hence $;(s2, yo) = 0 on the plasma 
boundary itself and also ¢,(0, yo) = 0 on the symmetry line (S5) from (5.2d). 

A pair of analogous expressions for integration in the y direction can also 


be written down for a given xo location, and these are 


y 
i Pyy (Xo, Y) dn = dy (xo, y) 


_ [" ( e@oy) . P¢ 
ct €0 ae: 


where s, is the height of the plasma boundary above the source cup at Zo, 


) dn, 0 S Lo < S1, 
XO,Y 


and where s; is the maximum extent in x that the boundary extends into the 
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region, and 
y 

He Pyy(Xo, Y) dn = Py(Xo, Y) 
0 


y , Oo? 
= 64(00,0)— f (saat ee 


€o 


dn, SX tg S Uy, 
L0,y 


(5.9a) 
where xp is the location of the right hand boundary, $3. It is noted that the 
term, ¢,(20,0) is present on the right hand side of (5.9a) because its value 
is not known there. 

A minor drawback with using these integral expressions (that was absent 
in the one-dimensional case) is the presence of the second derivative terms, 
byy(%, yo) and dz2(xo, y), on the right hand sides of (5.8) and (5.9), respec- 
tively. To successfully incorporate (5.7) and (5.9) into an iterative scheme, 
we must determine values for these terms. This was avoided in the one- 
dimensional formulation, but nevertheless does not present a problem in the 


strategy proposed here. 


The system of equations, (5.1), written in the required two-dimensional 


Cartesian form can now grouped together as 


OVz Oy L OVr = af (aa 4 CMY 
™ Ja 


iy — 


Ox ” Ox ot ay Oy? 


dé, (5.10a) 
&,YO 


du, Ovuy  Ovy gq / Yl p(xo,y) Pd do 
“y Oy "y Oy Oot ™m | 7 €9 = OFF a. MT OY || 
(5.10b) 
O(pvz) wa O(pry) _ _ Op 
An Bij BE and, (5.10c) 
2 2 
OO ODA, DP (5.10d) 


where a and 6 refer to the integration starting locations in x and y, respec- 
tively. As described above, these are either s,, or 0, or s, or 0 depending upon 
where the integral is to be determined from; the term ¢,(%o, a) (in (5.10b)) 
is only present if b = 0. It is noted that the divergence terms in (5.1) have 


been expanded (in a Cartesian manner) to give (5.10c) and (5.10d). 


5.1.1.1 Dimensionless Cartesian Form 


The system, (5.10), is now written in a dimensionless form, and to do this, 


as in §4.3.1, we define the dimensionless variables 


o= ~. (5.11a) 
 : 

— ee (5.11b) 
. Uz 

te= (5.11¢) 
& UV. 

i= es (5.11d) 

t= = (5.11e) 

72 a (5.11) 

and, f= (5.11g) 


Here, L is the largest physical dimension (% and y are equally scaled) and 
Om: Pm: Um and t,, are as defined in §4.3.1. In particular, we restate v,, and 


tim as 


(5.12) 


tm = —. (5.13) 
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The Dimensionless Gauss’ Law Since the integral expressions on the 
right hand side of (5.10a) and (5.10b) are determined from Gauss law, (5.10d), 
we write its dimensionless form first. The dimensionless variables particular 
to Gauss’ Law, as written above, are those used in §4.3.1, and so we simply 
state its dimensionless form as 


Po Ad 2 
DR? + ay? = —ay1/p, (5.14) 


where a, is given by (4.32) and similarly takes a typical value O(10). The 


integrals, (5.7) and (5.9) can therefore be written in dimensionless form as 


oo) sf? ff 2 OG 
We o = -{ ai p(z, Yo) + ag 7 dé, and 
z,Yo v,YO 
GON PP noe. oe OD de 
a| -~f onp(Gos9) + aaa] HS | <3 
LOY £0,Y £0,b 


where @ and 6 are the dimensionless integration starting locations; these are 


analogous to a and 6 in (5.10a) and (5.10b). 


The Dimensionless Lorentz Equations By substituting (5.15), along 


with the relevant dimensionless variables, into (5.10a) and (5.10b), we obtain 


. Wi,  . Wy Be DU eins cue, SOOO 
eae + Oyae + Br = af a1 p(X, Yo) a2 = ey dé, (5.162) 
#90 
and 
_ 0,  _ Ody Diy if ad 
ee Saye ao ; D~D d Ae 9 
U Oy Vy Oy 7 Ot a2 i: a1 p(Xo y) + Ox2 _ 1) Og me 
XLO,Y xO, 
(5.16b) 
where a2 is given by 
OL 
at 


from (5.12) and (5.13). The parameter a» has a typical value of ~ 107° and 
will tend to remove the effects of the integrals on the right hand sides of 
(5.16a) and (5.16b), unless the second derivative terms are significant. This 
will tend to be true just away from the plasma boundary where the curvature 


in @ is at its greatest. 


The Dimensionless Continuity Equation Substituting the relevant di- 


mensionless variables into (5.10c), we obtain 


Prim ( O(Ptx) | O(Pey)\ _ Pm OP 
L ( oa Og ) tm Ot’ on 
and by noting that t,, is given by (5.13), (5.17) can be written 
ONE) 50a) Oe (5.18) 


Dimensionless Boundary Conditions To be consistent with the dimen- 
sionless system, (5.14), (5.16) and (5.18), the boundary conditions, (5.2), 


must also be restated in dimensionless form and these are 


OE 
x)= 
= os, #>0,%€5,US, (5.19a) 
3 OT 
x,t) = — 
(x, t) bn 
= or, #>0,%€S, (5.19b) 
V 0(x, ft) - ay, = 0, i>0,x€S, (5.19c) 
dd : rane 
| =o i>0,x€ 5s (5.19d) 
x,t 
A(x) = mt i>0,x€5; (5.19e) 
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p=0, t>0, & € So, (5.19f) 


6-0, t=0,xER, X€¢S,USy, (5.19g) 
ap F oa 
| =0, #>0,xES 5.19h 
Ae x 5 ( ) 
s Vien 
v=—n, 
Um 
= 6;(%)h1, t>0,xES, (5.193) 
v=0, #=0,%ER, x¢S; (5.19}) 
v- iis = 0, bs Use Sa (5.19k) 


Here, R is the dimensionless region enclosed by the boundaries S, to Ss; 
the vectors fi, and fis are the unit normals to S; and $5, respectively. The 
dimensionless gradient operator, V, in (5.19c) has been written in Cartesian 
form in (5.19d) and (5.19h), where it is clear that potential or charge density 
do not change across the line of symmetry, at = 0. Typically, ds = 1 and 


dr = 0, since dr = 0 and dm = os. 


5.2 Strategy for the Determination of the Two- 
dimensional Steady-State Solution 


In a vacuum tube, in which the emission boundary is fixed (such as a vacuum 
tube diode), high emission current densities accompanied by low acceleration 
voltages will give rise to a charge build-up immediately above the emission 
surface. This occurs because the high charge flux from the emission surface 
is not matched by the acceleration of charge away from the surface, as a 


consequence of the relatively low acceleration potential. The effect of such 
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a build-up of space charge is that a steady condition ensues in which charge 
leaving the emission surface behind the space charge build-up is repelled by 
it ([3], p191). This condition is known as space charge limitation and can be 
modelled in one dimension by the Child-Langmuir relationship! ([3], p171). 

If the space charge limitation is severe, the electric field at the build-up 
of charge will be neutralised, effectively indicating the location of the new 
emission boundary ([3], p191). The numerical method presented in this chap- 
ter exploits this property to determine the location of the plasma boundary. 
It does this by calculating a fixed boundary solution, only terminating this 
calculation when either a maximum in potential exists within the solution 
region”, or when the numerical method begins to become unstable, as a con- 
sequence of the space charge limitation. Instability is identified by examining 
the mean of the relative root mean square change in solution, of all solution 
variables (called MRMS in ensuing paragraphs), at each time-step. Generally 
convergence is observed in that there is a rapid decrease in MRMS with time- 
step, until space charge limitation becomes apparent when MRMS generally 
begins to increase. This is a phenomenon also observed when using other, 
off-the-shelf, fixed boundary ion beam codes (for example, OPERA2d), that 
calculate the fixed boundary solution in a slightly different way by tracking 
particles individually in a Lagrangian manner. When space charge limitation 
occurs with these codes, fixed boundary solution convergence is usually not 


'This relationship is effectively the fixed boundary equivalent solution to the planar 


free-boundary solution developed in Chapter 3. 
? At this point the charge entering the region is not being accelerated away fast enough 


by the accelerating potential, resulting in space charge build-up and current limitation, 


and a consequent maximum in electrostatic potential ([3]). 
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achieved. 

When either a potential maximum within the solution region or an in- 
crease in MRMS (or both) has been identified, the solution is examined and 
the new boundary location determined to be the location of the zero electric 
field level curve (satisfying (5.19c)) near the existing fixed boundary. Further 
iterations can then be performed, with the boundary location generated pre- 
viously forming the new emission surface. This process can be repeated until 
either the boundary does not move any further, or until there is no further 
build-up of space charge immediately above the emission boundary; both 
conditions should occur at the same time. At this point, the free-boundary 
location has been determined. 

This solution strategy has been implemented in Matlab (as with all pre- 
vious models), and we now describe the solution algorithm, as implemented, 


in more detail. 


5.2.1 Solution Algorithm for the Two-dimensional Steady- 
state Calculation 

The method proposed in §5.2, as implemented, can be broken down into the 

following procedure. The algorithm steps are as follows. 

Step 1 


e Set up the fixed solution region with n, equally spaced nodes in the x 
direction, and n, equally spaced node in the y directions; the spacing 
in each direction is labelled AZ and Ay, respectively, where Az = 


Ep/(nz — 1) and Ay = gr/(ny, — 1). As in the one-dimensional time- 
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dependent model, the subscripts 7 and j are used to represent i*® node 
in Z and the j" node in g. Ensure the right hand edge (Zz) of the fixed 


emission region coincides with a nodal point. 


e Set the time-step size At. As in the one-dimensional time-dependent 
model, the superscript k is used to represent a variable at the k‘" time- 


step. 


e The fixed emission boundary sits adjacent to the line of symmetry (S5 
in Figure 5.1), along the % axis; it is labelled $; = Sf(%,0) here, with 


its right hand edge being located at (%,, 0). 


Step 2 


e Initialise the solution variable matrices for the variables d, Pp, Vz, and 
vy. The solution is stored at each nodal point in matrices of dimension 


(ny,Nz), accounting for each nodal point within the region. 


e Add any constant (in time) boundary conditions to these variables, 


prior to beginning calculations, as follows. 


— The charge density (at each time-step) is ramped, in space, along 
the emission boundary to avoid a discontinuity at the right hand 


edge; it takes the form 


A(é) = " + exp (oe =) yo 


where here, c,/L is the midpoint of the ramp transition and w,/L 
is a measure of the ramp transition width (previously described 
in §4.3.2.3). This gives rise to a curve similar to that shown in 


Figure 5.2. 
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Figure 5.2: Constant charge density profile along the emission boundary. The 
region line of symmetry is located at x/L = 0, the midpoint of the ramp transition 
is at c,/L (in this case at x/L = 0.16, shown by the black square marker) and a 
measure of the transition width (the distance between the two red, circular markers) 


is given by w;/L (w,/L = 0.0286 in this case). 


— The velocity magnitude (at each time-step) along the boundary, 
0;(Z), has a similar ramp applied. In the fixed boundary case, 
since ions are emitted normally away from the emission boundary, 


only the component 0, is affected. 


— The conditions, (5.19a) and (5.19b), are applied to the potential 
solution matrix. 
Step 3 
e As a starting point, calculate the equivalent one-dimensional steady- 
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state potential solution for the boundary conditions applied, using the 
analytic expression, (3.28); map this solution onto the line of symmetry 
in both the potential and charge density solution variables (the charge 


density being calculated from (4.1c). 


e Using the one-dimensional solution along this boundary, in conjunction 
with the other boundary conditions already applied, perform relaxation 
on the potential and charge density solution variables. The relaxation 
process applied is effectively the same as seeking the solution of the 


Laplace equations, 


Ld &d 

OE am Ba? = 0 and 
OD, Op _ 

O52" Oy” 


for the potential and charge density, respectively. It is done using a 
five point nodal averaging procedure, where for a general point, bij; 


say, within the region, the (k + 1)" average would be given by 


PEP = (big + Chay + Cijar + Gj) /4- 


Since we are only seeking a steady-state solution, the initial solutions in db 
and p, generated by the relaxation procedure, are a suitable starting location 


for the steady state calculation; 
Begin iterations on the fixed region, proceeding in the following order. 


Step 4 
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e Calculate baa from the most recent potential solution, using the second 


order central difference, (4.55). Calculate the updated 0, component, 
RAY 5 ke x ms gS on aN 
Vel; j a Vel; 5 + At 4 a2 OP; gr bya}. | dé 
(0) ay) 


_ O,]* [_ d6,]* 
~ Se] - paz]. (5.21) 
J 


from (5.16a), using upwinded differences (as in the one-dimensional 


case) for the first order derivatives. The integral in (5.21) is calculated 
using the trapezium rule. 
Step 5 


e Calculate ¢;; from the most recent potential solution, using the second 


order central difference, (4.55). Calculate the updated 0, component, 
~ (k+l ~ 1k 4 v xk rah 
dyl;,; = dy |; 5 + At {a2 | (art, + bz% ) dé 
0 uJ 


_ Ob,]**" [_ d,]* 


a,j a,j 


from (5.16b), using upwinded differences (as in the one-dimensional 
case) for the first order derivatives, and using the 0, component most re- 
cently calculated. The integral in (5.22) is calculated using the trapez- 


ium rule. 


Step 6 


e Calculate the updated charge density, 
oe is, O 
~ktl] __ xk SO (xk = jk+1 Oo (xk ~ k+l 
Pig = Pi3 —At x (a, Vel; j ) a ag Ga dyl 5; )} (5.23) 
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from (5.18), using the most recently calculated values for 0, and dy. 
Derivatives in (5.23) are calculated using first order upwinded differ- 


ences, as in the one-dimensional model. 


Step 7 


e Solve Gauss’ Law, (5.14), using the most recently calculated values for 
p. This is done using a standard, second order five point difference 
stencil (as described in [26], Chp 6). Examine the resulting solution 
and if a potential maximum exists within the solution region (instead 
of at the plasma boundary), terminate the fixed-boundary calculation 
at the current time-step. The Matlab “find” command is a suitable tool 


for performing the examination. 


Step 8 


e Generate the MRMS at the current time-step. This is done by calcu- 
lating the RMS change in solution (between the current and previous 
time-step) as a fraction of the RMS solution at the current time-step, 
for each of the four solution variables, v,;, vy, p, and é. For p, for 


example, the relative RMS change in solution is calculated as 


“sh As 2 
a Pap ey) 
a) 


RMS(p) = SA —__., 
"TE Seuny 
ig 


where the product (n,n,) is the total number of nodes. Once this is 
known for all solution variables, the mean of all four (previously called 


the MRMS) is calculated. 
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Step 9 


e Examine the rate of change of MRMS and if this has a positive mean 
over, say, the previous ten time-steps, terminate the fixed boundary cal- 
culation at the current time-step. If this condition is not true, continue 


the iterations by returning to Step 4, above. 


Step 10 


e Once the initial fixed boundary calculation has been concluded, ex- 


amine the potential solution and calculate the electric field magnitude 


a-@ 


e Determine the level curve corresponding to an electric field magnitude 


from 


of zero. This is done using the Matlab “contourc” command, with a 


requirement for a single contour at height zero. 


e Identify those fixed nodes, within the solution domain, that are clos- 
est to the zero electric field magnitude level curve. These nodes form 
the location of the plasma boundary for the next stage of calculation. 
An example of the result of this process is shown in Figure 5.3. The 
most convenient tool for locating the nearest nodes is the Matlab “find” 


command. 


Step 11 


e Calculate the gradient of the new boundary at each node along it; from 


the gradient, calculate angle of the inward (into the region) directed 


169 


1 lo) o o o o o 
0.9% ° ° ° ° fe) ° fo} ° ° ° 4 
° ° ° ° fo) ° ° ° fo) fo) 
0.87 
° ° ° ° fe} ° ° ° ° ° 
0.7+ 
° fo) ° ° ° ° ° fo) ° oO 
© 0.6- 
£ 
5 ° fe} fe} fe} te) fe) ° ° fe} fe} 
5 0.57 
° ° ° ° ° ° ° ° ° ° fo) 
a 
> 04+ 
° ° fo) ° ° ° ° ° fo) fo) 


° ° ° fe} 
oO ° ° ° 
oO ° fe} fe} 


! ! ! ! ! ! ! ! 
0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 
X location (m) 


Figure 5.3: The location of the level curve corresponding to zero electric field 
magnitude, and nearest grid nodes. The blue line is the level curve, whilst the black, 
square markers are the nodes identified as being nearest this curve. Information 
is displayed in dimensional space, with the right hand edge of the flat emission 
boundary (ag) being located at x = 0.01m. After the initial fixed boundary 


solution, nodes below the blue line are considered to be outside of the region. 


normal at those nodes. If the angle of the inward normal at the p'* 
node on the new plasma boundary is labelled 6,(X,), then the velocity 
boundary condition, (5.19i), as applied to the 0, and d, velocity com- 
ponent solution variables along this boundary, is given by the pair of 


expressions 


Ux(X») = 0; cos(6,(X,)), and 


e Apply the charge density condition, (5.19e), along the new boundary. 


Step 12 


e Repeat the iterative process described above, but when necessary, in- 
tegrate in the % and ¥ directions from the expanding boundary (not 
including the nodes behind it); remove the rows and columns, corre- 
sponding to these nodes, from the matrix generated within the Gauss 
Law Poisson solver, moving those values on the new boundary to the 
source vector on the right hand side of the system of equations gen- 
erated by the Poisson solver. Calculate MRMS at each time-step of 
the iterative procedure. Halt the procedure if a maximum in potential 
exists inside the region, or if MRMS begins to increase; re-examine the 


solution. 


5.2.2 Test Case 


To test the method, a set of conditions that will cause plasma to expand 
into the solution region, but not cause the emission region to become too 
distended, are applied. Generally, the application of a relatively high current 
density at the emission boundary for a normal operating voltage will cause 
this to happen. Indeed, usual operating conditions, applied to a neutron 
tube, should cause some bulging of the boundary into the solution region, 
and so we choose to apply a source potential of 120kV, a source current 
density of 18800Am~°, and an initial particle energy of 50eV. The applied 
conditions detailed here are listed in Table 5.1, along with other numerical 


parameters. 
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56 Number of nodes in x direction. 
1/1100 m_ | Nodal spacing in x direction. 


112 Number of nodes in y direction. 


1/11100 m_ |} Nodal spacing in y direction. 
1x 107-!'s | Time step size. 
0.05 m Region width (length scaling parameter). 
0.01 m Nominal source-target spacing. 
0.01 m Width of emission region. 
100 kV Source potential. 


OkV Target potential. 


18800 Am~? | Source current density. 


50 eV Initial particle energy. 


Table 5.1: Numerical parameters for the two-dimensional steady-state test prob- 


lem. 


5.2.2.1 Results 


Selected results are displayed for three stages of calculation, in which the 


boundary is seen to advance from the initial fixed location. 


Fixed boundary solution - stage 1 After exiting the initial calculation 
stage, the calculated boundary location was found and this is shown in Fig- 
ure 5.4. Here, because the new calculated boundary location is so close to 


the initial, fixed emission region, the nodes within the region that are closest 
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to the new boundary, are still located along the initial fixed boundary. The 
calculated new boundary location is shown as the blue line, and the closest 


nodes as the black, square markers. Figure 5.5 shows the modified boundary 
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Figure 5.4: After the initial calculation stage, the new boundary location is cal- 
culated to be the blue line. The nodes within the fixed grid that are closest to this 
boundary are shown as the black, square markers. A small part of the underlying 


computational grid is also shown. 


conditions applied to the v, and v, variables, respectively, as a result of the 
shape of the new boundary shown in Figure 5.4. Whilst, on the computa- 
tional grid, the boundary isn’t changed for the next stage of the calculation 
(since the nearest nodes are those on the existing fixed boundary), the v, 
and vy components, along the boundary, are updated. It can be seen that 
the y directed velocity component along the boundary is some one hundred 
times the magnitude of the x directed component. This (the v, boundary 
condition) will have very little effect on the next stage of the calculation. 


The potential solution at this stage is shown in Figures 5.6 and 5.7, where 
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Figure 5.5: The modified boundary conditions for vz and vy as a consequence of 


the new boundary shown in Figure 5.4. 


in Figure 5.6, the build-up of charge ahead of the fixed emission boundary 
is causing, what is now, typical curvature in the potential solution there. 
Figure 5.7 shows potential level curves (40 curves spaced 3kV apart), where 
a "bulge" immediately above the emission region is apparent. 

Figures 5.8 and 5.9 show the charge density solution at this stage. In 
Figure 5.8, the ramped, constant charge density across the emission region 
can be seen; it changes from a maximum of jp/v; = 0.2717 Cm~® to zero 
rapidly at 4/5 distance along (traversing from right to left in the upper right 
hand part of the figure). Figure 5.9 shows the corresponding charge density 
level curves (40 curves, equally spaced from 0 to 0.2717 Cm~%). 

Figures 5.10 and 5.11 show the x component of velocity, v,. The sharp 


ridge in Figure 5.10 indicates that particles, in this location, are being accel- 
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Figure 5.6: The initial potential solution surface. 
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Figure 5.7: Initial potential solution level curves. The calculated, new plasma 


boundary is just visible and is shown as the blue line to the bottom left of the 


figure. 


erated outwards as they travel towards the target. Physically, this is caused 


by mutual charge repulsion within the ion beam. A slight, overall beam 
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Figure 5.8: The initial charge density solution surface. The ramp along the emis- 
sion surface (to the upper right of the figure) is visible, where the charge density 
at the emission surface decays rapidly from its maximum (jo/v;) to zero when at 


4/5 distance along. 
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Figure 5.9: Initial charge density level curves. The highest charge density is shown 


in red, immediately above the emission region. 
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Figure 5.10: x component of velocity. The sharp ridge is evidence of mutual 
charge repulsion within the ion beam and shows evidence of lateral acceleration 


caused by this. 


splaying is visible in Figure 5.11; at this accelerating voltage, experimental 


experience indicates that only a small degree of beam splaying would occur. 


For completeness, Figure 5.12 shows the y component of velocity. The 
curved shape of v, as a function of source-target distance was seen in the one- 
dimensional model; for comparison, the expected one-dimensional velocity 
of a deuteron under acceleration from the steady-state potential solution is 
shown in Figure 5.13, along with the two-dimensional v, component, along 
the y-axis (or line of symmetry). In both cases, the same ultimate speed at 
the target is reached (indicating some consistency in calculation), but the 


two-dimensional curve is distinctly more curved. 
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Figure 5.11: x component of velocity level curves. A small degree of beam splaying 


is evident. 
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Figure 5.12: y component of velocity. 


Figure 5.14 shows the MRMS (described above) for the 400 iterations of 
this stage of calculation. There is a clear reduction in MRMS as the calcu- 
lation proceeds, indicating that until the stage was terminated, the solution 


was converging. 
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Figure 5.13: Comparison between two-dimensional vy component (along the y- 
axis line of symmetry) and the one-dimensional calculated velocity of a deuteron 
under the acceleration caused by the calculated steady-state accelerating voltage. 
It is noted that the one-dimensional curve does not reach the source boundary, as 


a consequence of the steady-state reduced region size. 
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Figure 5.14: MRMS with iteration number - stage 1. 
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Modified boundary solution - stage 2. The next fixed boundary calcu- 
lation stage proceeded for a further 375 iterations before termination. Whilst 
the boundary calculated in stage 1 did not result in any change in the region 
size, velocity components were modified as mentioned above. This, second 
stage, resulted in more significant change in boundary location as shown in 
Figure 5.15. Here, the new boundary is shown as the blue line, with the 
underlying region nodes that are nearest, shown as the black, square mark- 
ers. For comparison, the red line indicates the location of the boundary from 


stage 1, above. Figure 5.16 shows the modified boundary conditions applied 
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Figure 5.15: Calculated boundary location for stage 2. The blue line shows the 
new boundary location and the black, square markers are the nodes nearest the 
boundary on the underlying grid. The red line is the boundary location calculated 


from the previous stage. 


to the uv, and v, variables, respectively, as a result of the shape of the new 


180 


boundary shown in Figure 5.15. It can be seen that the y directed velocity 
component along the boundary is some one hundred times the magnitude 
of the x directed component. The peak x directed velocity component at 
the boundary is now some twenty times smaller than the y directed velocity 


component. 
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Figure 5.16: The modified boundary conditions for v; and vy as a consequence of 


the new boundary shown in Figure 5.15. 


Figures 5.17 and 5.18 show the potential solution at stage 2. Within 
Figure 5.17, nodes behind the new boundary have been removed to highlight 
the new emission boundary (the missing nodes can just be seen at the upper 
right of the figure) and will not contribute to the following calculation stage. 
The calculated emission boundary is shown in Figure 5.18 and whilst the 
potential solution does not look significantly different to that shown in stage 


1 (Figures 5.6 and 5.7), examination of the level curves (40 curves spaced 
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3 kV apart in both cases) near the emission boundary in both cases show that 
there is some boundary advancement (i.e. the two curves nearest the emission 
region are advancing from stage 1 to stage 2). This indicates advancement 


of the plasma boundary. 


New boundary 
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x location (m) 
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Figure 5.17: Stage 2 potential solution surface. Nodes behind the new emission 
boundary have been removed to highlight their location (just visible to the upper 


right) and do not contribute to the next stage. 


As a way of determining the extent of the beam splaying at this stage, the 
charge density solution was scaled (by the reciprocal of the particle charge) 
in order to determine the particle number density across the region. Figure 
5.19 shows the particle number density across the target (a logarithmic scale 
is used on the ordinate axis), where it can be seen that after ~ 0.015 m, 
there is less than one particle per unit metre within the region. Therefore 
we conclude that the very outer edge of the ion beam strikes the target at 
the location « = 0.015 m. This indicates some degree of beam splaying, 


considering that the majority of charge is emitted from a region between the 
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Figure 5.18: Stage 2 potential solution level curves. There is advancement of the 
curves nearest the boundary from stage 1. The new, calculated plasma boundary 


is shown as the blue line to the bottom left of the figure. 


y axis and x = 0.008 m. 

Figure 5.20 shows MRMS for the combined stages one and two (775 iter- 
ations). The discontinuity at iteration number 400 is likely to be caused by 
the minor modifications introduced as a consequence of the boundary shape 
in stage 1 (i.e. the velocity component modifications). Again, there is a clear 
reduction in MRMS as the calculation proceeds, indicating that until the 


stage was terminated, the solution was converging. 


Modified boundary solution - stage 3. The final stage of calculation 
proceeded for a further 350 iterations, where it was terminated. Figure 5.21 
shows the boundary location at the end of this stage. As before, the new 
boundary is shown as the blue line, with the underlying region nodes that 


are nearest, shown as the black, square markers. For comparison, the green 
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Figure 5.19: The particle number density across the target. The very outer edge 
of the beam is considered to be when the number density falls below one particle 


per cubic metre (in this case at 0.015 m). 
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Figure 5.20: MRMS with iteration number - stages 1 and 2. 


and red lines indicate the location of the boundary from stages 1 and 2, re- 


spectively. For comparison, the location of the boundary calculated from the 
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one-dimensional planar, analytic solution ((3.23)), for the given set of con- 


ditions is shown. It is noted that part of the boundary (x > 0.0063 m) does 
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Figure 5.21: Calculated boundary location for stage 3. As before, the blue line 
shows the new boundary location and the black, square markers are the nodes 
nearest the boundary on the underlying grid. The red and green lines are the 


boundary locations from steps 2 and 1, respectively. 


not penetrate into the body of the region. This is due to the applied source 
charge density ramp decaying rapidly beyond this point. In this region, the 
charge density is relatively low, and combined with the accelerating voltage, 
will give rise to a recessed boundary. Figures 5.22 and 5.23 show the po- 
tential solution at stage 3. As before, within Figure 5.17, nodes behind the 
new boundary have been removed to highlight the new emission boundary. 
The calculated emission boundary is shown in Figure 5.23, where the blue 
line to the bottom left hand corner indicates the plasma boundary location. 


The gap between the boundary and the first level curve is an indication that 
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Figure 5.22: Stage 3 potential solution surface. Nodes behind the new emission 


boundary have been removed to highlight their location. 


the potential solution is not changing rapidly there ((5.2c)). Indeed Figure 
5.24 shows the calculated electric field magnitude across the region, where it 
can be clearly seen that the field magnitude is zero at the boundary surface 
(shown as the blue line). 

Figure 5.25 shows the particle number density across the target. Here, 
the very outer edge of the beam, based on the previous definition, is found to 
be at about 0.014 m, indicating that the beam has focussed inwards slightly 
from the previous stage, this is slightly counter-intuitive, but the previous 
stage was an interim stage of the calculation. 

Finally, figure 5.26 shows MRMS for the combined stages one, two and 
three (1125 iterations). The quite large discontinuity at iteration number 
775 is most likely to have been caused by the boundary modifications intro- 


duced as a consequence of the boundary shape in stage 2 (i.e. the velocity 
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Figure 5.23: Stage 2 potential solution level curves. The calculated plasma bound- 
ary is shown as the blue line to the bottom left of the figure. The gap between the 


boundary and the first curve indicates that the solution is not changing rapidly 


there. 


component modifications and the movement of the charge density conditions 
to the distended boundary). However, once again there is a clear reduction 
in MRMS as the calculation proceeds, indicating that until the calculation 


was terminated, the solution was converging. 


5.3. Summary and Conclusion to the Two-dimensional 


Steady State Model 


A strategy for a two-dimensional, planar steady-state solution to the plasma 
boundary problem has been proposed, and this is based in part on the sec- 


ond one-dimensional, time-dependent numerical method. The method uses 
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Figure 5.24: The electric field magnitude across the region. The electric field 
reaches zero at the plasma boundary, as required by (5.2c) (or the dimensionless 


equivalent). 


the notion of space-charge limitation, whereby charge builds up away from 
a fixed boundary on the interior of an acceleration region, causing a poten- 
tial maximum and drop in electric field. This happens as a consequence of 
the fixed boundary being unable to move, and thus for a given region, ap- 
plied electric field and source charge density combination, charge cannot be 
accelerated away from the fixed boundary sufficiently quickly. In a moving 
boundary case, this would result in the boundary moving to increase the 
electric field near its surface, thus allowing charge to be accelerated away. In 
the fixed boundary case, however, the build-up of charge effectively creates 
a new boundary, whereby the electric field at its surface is zero. We detect 
the location of this boundary, by looking for a potential maximum within 


the region. 
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Figure 5.25: The particle number density across the target. The very outer edge 
of the beam is considered to be when the number density falls below one particle 


per cubic metre (in this case at 0.014 m). 
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Figure 5.26: MRMS with iteration number - stages 1, 2 and 3. 


Within the numerical method proposed, this is done in a number of dis- 
crete stages (three in the test problem shown here), where charge density 


and velocity boundary conditions are moved from the previous boundary lo- 
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cation, to the new boundary location, and iterations continued until a further 
potential maximum is found within the region. This process continues until 
the measure of solution variable change (MRMS) is deemed to be changing 
sufficiently slowly that the solution has converged. 

The test case here shows that the method is entirely plausible and whilst 
not exhaustive, does indeed show the method converging (in three separate 
stages), whereby the final resting place of the boundary is close to that of 
the one-dimensional equivalent solution. The method, as implemented, does 


not allow for recessed boundaries, and dealing with this is future work. 
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Chapter 6 


Summary, Conclusions and 


Further Work 


6.1 Summary and Conclusions 


The problem of determining the plasma free boundary location within an ion 
accelerator set-up has been considered for many different applications, and 
is detailed in 81.1.2. The work presented in this thesis has concentrated not 
only on determining the solution to the free boundary problem originating 
in a neutron tube acceleration region, but also on the full time-dependent 
moving boundary problem, of which the free-boundary problem is the steady 
case. 

The formulation of system of non-linear equations describing the moving 
boundary problem is detailed in Chapter 2, with the full system of non-linear 
equations describing the time-dependent problem being given by (2.13), and 
the reduced system of non-linear equations describing the steady-state prob- 


lem being given by (2.19). 
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6.1.1 One-Dimensional Steady-state Problem 


The one-dimensional free boundary problem is considered in Chapter 3, 
where the full steady-state system, (2.19), is reduced to one dimension, giv- 
ing rise to the planar system, (3.1) (with the associated boundary conditions 
(3.3)), and the radial system, (3.30) (with the associated boundary condi- 
tions (3.29)). The full analytical one-dimensional solution to (3.1), subject 
to (3.3) is developed and is eventually given by the rather unwieldy equation 
(3.28). However, the similar radial problem, detailed in §3.2, does not appear 
to be analytically soluble, and so a numerical method is developed that is 
tested against the analytical planar solution. This is successfully extended to 
the radial problem, where solution features that are particular to the radial 


problem, are highlighted. 


6.1.2 One-Dimensional Time-dependent Solution 


In Chapter 4, a one-dimensional time-dependent solution to the the full sys- 
tem of equations, (2.13), is sought. This is achieved by initially reducing 
the complexity of the system by making the physically realistic assumption 
that any magnetic field influence (either by that applied externally, or by 
that induced as a consequence of the charged particle flow) can be neglected. 
By doing this, the full one-dimensional time-dependent system can then be 


written as (4.1), with the associated boundary conditions, (4.2). 


The one-dimensional time-dependent system is analysed, giving rise to the 
analytic expressions, (4.9), describing characteristic curves in the charge free 


region ahead of the advancing wavefront. Additionally, manipulation of the 
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time-dependent system allows it to be reduced to the parametric, ordinary, 
Riccati differential equation, (4.18), in terms of the solution variable Ov/Ox 
(the rate of change of velocity with distance), as a function of the parametric 
coordinate o measured along characteristic curves defined within the solu- 
tion region. The form of the solution to this equation can be relatively easily 
determined to be (4.23), and is written in terms of Ov/Ox as a function of 
time (to within two constants). One of the resulting constants can be rel- 
atively easily determined, and it is also believed that the other can also be 
determined, although the resulting equation for 0v/0z is likely to be implicit. 


This is something that will be considered in due course. 


In order to develop a full time-dependent solution to the one-dimensional 
time-dependent system, (4.1), it is first re-written as the dimensionless equa- 
tions (4.31), (4.33) and (4.34), with the associated dimensionless conditions, 


(4.36). Following this, two approaches to solving the system are proposed. 


6.1.2.1 Numerical Method 1 for the One-dimensional Time-dependent 


System 


The first method proposed for solving the system, (4.31), (4.33) and (4.34) 
(with associated boundary conditions), begins by mapping the system of 
equations that are to be solved on a moving domain, to a logical, fixed do- 
main, where standard finite differences can be used to represent the deriva- 
tives within the now, logical, system. The mapping process incorporates a 
solution dependent spatial mapping function which controls nodal spacing in 


the physical region, and which also varies according to the changing physical 
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region size. Additionally, the physical time-step is also adjusted by a solution 
dependent time-step controlling function function, which allows non-constant 
physical time stepping to be achieved. Once the time-dependent solution has 
been determined at equally spaced nodes and at equal time-steps in the log- 
ical region, it can be transposed onto the unequally spaced nodes and time- 
steps in the dimensionless physical region, with the full dimensional solution 
being determined accordingly. Included within the method is a mechanism 
for determining the location of the moving boundary in dimensionless phys- 
ical space. 

A numerical algorithm for determining the solution of the one-dimensional 
time-dependent problem in the logical domain is developed and listed in 
detail; this has been implemented using Matlab. Numerical experimen- 
tation indicated that whilst field solutions generated by this method look 
viable, in order to calculate the boundary location using the mechanism 
proposed, a small, multiplicative parameter must be introduced to reduce 
calculated boundary movement at each time-step; without such a parameter, 
the method rapidly fails. The use of this parameter renders the calculated, 
time-dependent boundary movement unphysical, whilst additionally causing 
the algorithm to execute very slowly. These things aside, calculation results 
are presented for a particular test case in which the charge density applied 
at the emission boundary is ramped with time!. Results show viable field 
solutions, and indicate that the solution dependent spatial and temporal 


mapping functions are operating as expected. However, as mentioned, the 


TAs aresult of the slow execution times, to examine the effects of the source charge den- 
sity switch-on, a rapid charge density ramp is applied; this is only just within experimental 


range, and is not typical of normal tube operation. 
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rate of boundary movement is entirely unphysical. 


6.1.2.2 Numerical Method 2 for the One-dimensional Time-dependent 


System 


To avoid the problems experienced in the first numerical method, a new ap- 
proach to the problem is proposed and developed. This new method reduces 
the number of equations in the system by one, by integrating Gauss’ law, 
(4.35a), and then incorporating it into the Lorentz equation, (4.35b), yield- 
ing the coupled pair of equations, (4.75). Additionally, further manipulation 
of Gauss’ law allows the integral equation, (4.79), to be written for the rate 
of change of boundary location. 

A numerical algorithm to implement the new approach is developed and 
listed in detail. Owing to much improved execution speed, a test case, in 
which a more experimentally realistic source charge density ramp, is exam- 
ined. Results clearly show the solution converging to a steady state, indi- 
cated by the fact that the calculated boundary location asymptotes to the 
expected analytic boundary location. Furthermore, boundary perturbation 
as a function of time, rapidly decays to zero as the calculation progresses. 
Comparisons between the exact, analytic solution and the time-dependent so- 
lution for differing numbers of nodes, show that as the nodal spacing across 
the region decreases within the time-dependent solution, the error in solution 


also decreases. 
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6.1.3. Two-Dimensional Steady-state Problem 


A strategy for solving the two-dimensional steady-state problem is proposed, 
in which space-charge limitation, caused a high source current density at a 
fixed emitter surface, causes the electric field at the location of the build-up 
of charge to be neutralised; this effectively indicates the location of the new 
plasma boundary. 

The strategy proposed uses an integral representation of the electric 
field (in a manner similar to that used in the second one-dimensional time- 
dependent method) within the Lorentz equation, (5.la). This equation, rep- 
resenting the time dependence of the velocity field, is then split into its two 
Cartesian vectorial components, yielding a pair of integro-differential equa- 
tions for the x and y velocity components. These equations, in addition to 
the Cartesian forms of the charge continuity equation, (5.10c), and Gauss’ 
Law, (5.10d), now represent the system to be solved for the two-dimensional 
steady-state problem. But as in the one-dimensional time-dependent model, 
they are initially written in dimensionless form. 

The numerical method developed, solves the system, (5.14), (5.16a), (5.16b), 
and (5.18), subject to the conditions, (5.19), on a fixed computational do- 
main, using differences to represent derivatives. The initial solution, with a 
flat, fixed emission boundary proceeds until a potential maximum is detected 
within the interior of the solution domain. At this point, iterations are ter- 
minated and the location of the zero electric field level curve determined; this 
location (or the computational nodes nearest it) now represents the emission 
boundary for the next stage of calculation. Velocity components along the 


boundary are then updated as a consequence of its, now, curved shape. The 


196 


charge density and potential boundary conditions are then applied to the 
new boundary and iterations restarted, with any integration from the emis- 
sion region starting at the new boundary; nodes behind the boundary are 
removed from the calculation. This process is repeated a number of times, 
until no further solution maximum is detected within the interior of the re- 
gion, and the rate of change of MRMS (the mean root-mean-square change 
in solution between iterations, of all four solution variables) falls below a 
specified tolerance. 

A test case, using typical conditions applied to a neutron tube (that it 
is expected will not result in a largely distended boundary), is examined. 
Results from the three stages of the calculation are displayed, where the 
boundary is clearly seen to advance into the solution region as iterations and 
stages progress. In stage 1, results show the advancement of charge into 
the region, with particles attaining a lateral acceleration as they progress 
towards the target. The MRMS at stage 1 clearly reduces with iteration 
number, showing convergence of the method until it was terminated. Results 
from the latter two calculation stages show the advancement of the boundary 
towards its resting place, just below the expected one-dimensional boundary 
location. An indication of the calculated beam splaying at the target is also 
shown, whilst MRMS decreases within each of the calculation stages. 

The method, whilst a little unrefined at this point, is clearly plausible, 
offering a viable solution approach to the steady-state two-dimensional prob- 
lem. Refinement of the method could include a better representation of the 
boundary, instead of it being represented by the discrete nature of the un- 


derlying computational grid. This could be done using a different numerical 
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approach (finite elements), or by modifying the difference equations near the 
boundary to accommodate its curved nature. Furthermore, as it stands, the 
method is unable to calculate recessed boundaries; part of the boundary cal- 


culated within the test case is likely to be recessed. 


The methods and results detailed in this thesis offer not only an analytic ex- 
pression for the planar one-dimensional steady-state plasma boundary prob- 
lem for particles emitted with an initial energy from the plasma boundary, 
but also an insight into the rapid, time-dependent charge flow and plasma 
boundary movement that occurs over short periods of time immediately fol- 
lowing neutron tube switch-on. These methods could be developed to model 


other charged particle accelerators that incorporate similar plasma sources. 


6.2. Further Work 


Throughout this work, a number of areas have come to light that merit 


further study; these are listed in their order of occurrence within the thesis. 


1. A full study into analytical solution to the one-dimensional time-dependent 
planar system of equations (presented in §4.2.2) is to be carried out, 
with the aim of determining closed form analytic solutions on charac- 


teristic curves throughout the charge infused solution region. 


2. A more accurate boundary representation within the two-dimensional 
steady-state model needs to be developed. This should accommodate 
situations in which recessed boundaries occur. 


3. Two-dimensional solutions are to be tested against an experimental 
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set-up designed to give real physical information on the dimensions of 
the charged ion beam, as it strikes the ion accelerator target. The ex- 
perimental set-up to be used is shown in Figure 6.1, where the source 
shield at the bottom left of the Figure is designed to represent planar 
emission with the incorporation of an elongated emission slot (this is 
labelled "Plasma Emission Region" in the accompanying schematic). 


The set-up shown is encased in a large (relative to the acceleration re- 
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Figure 6.1: Experimental accelerator set-up. The components shown in the 
schematic are encased in a sealed envelope, with the required potential difference 
between the source shield and target scintillator applied with a large pulse forming 


network. 


gion shown) sealed envelope, with the required short pulse acceleration 
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potential difference, applied across the acceleration gap between the 
source shield and scintillator target, by a large pulse forming network. 
Ions released from the plasma emission region strike the scintillator tar- 
get where they are absorbed, and their kinetic energy released in the 
form of light that can be photographed by a high speed camera. Com- 
parisons can then be made between the predicted ion beam width as 
it strikes the target, and the recorded scintillation photographs. Com- 
ponents for this experiment have already been designed and procured 


for use in existing facilities. 


This experimental study will determine the validity of the mathemati- 
cal model used and perhaps indicate modifications that may be neces- 
sary to capture the physics more effectively (such as tuning the charge 
density ramp in both space and time at the plasma emission region). 
Similar experimental studies, representing a cylindrically symmetric sit- 
uation, say, will also be helpful in the longer term in further validating 


the modelling methodology when used in associated problems. 


. The work presented here is potentially informative in the modelling of 
other related problems, and the methodology for solving them. Indus- 
try largely uses a Lagrangian approach to charged particle beam mod- 
elling, but we have shown that the Eulerian framework is equally viable, 
offering potential advantages not seen in the Lagrangian approach (for 
example short timescale resolution of boundary movement). So there 
is the prospect that a significant amount of research could follow from 


the start made here. 
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Appendix A 


D-T Fusion Reaction 


Cross-Sections 


For a fast moving particle incident upon a stationary target particle, the 
total fusion reaction cross-section is the apparent target area, as seen by the 
incident particle as it approaches the other. It is effectively a measure of 
the probability that a reaction will take place between those particles and is 
dependent upon the kinetic energy of the interaction. 

For a fusion reaction between two ions, if the incident ion strikes the tar- 
get ion with a low kinetic energy, it will be unable to penetrate the coulom- 
bic barrier existing between the two particle nuclei, thus causing the reaction 
probability (or the effective target area) to be low!. If the incident ion strikes 
the target ion with a very high kinetic energy, whilst it will be able to pene- 
trate the coulombic barrier existing between their two nuclei, the incident ion 


Tt is understood that the larger the target ion apparently is, the more likely the incident 


ion is to hit that target. 
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will be less likely to be captured by the target ion nucleus’, thus scattering 
from it and again causing the reaction probability to be low. At a range 
of energies specific to the given reacting ion species, the probability of the 
incident ion nucleus being captured by the target ion nucleus is relatively 
high; at these energies, fusion between the particles is more likely. 

We now explain how, for a fast moving ion incident upon a target ion 
that is stationary relative to the laboratory, the reaction cross-section in the 
laboratory can be determined from that given in, what is known as, the 
centre-of-mass (CM) frame of reference. The following analysis is taken, in 


part, from [21]. 


A.1 Kinetic Energy in the CM Frame of Ref- 
erence. 


For a pair of reacting particles, the laboratory frame of reference (L) is that 
frame which is stationary with respect to the laboratory; the CM frame of 
reference is defined as the frame in which the total momentum of the reacting 
particles is zero. To determine the relationship between the total kinetic 
energy of such a pair of particles measured in the CM frame of reference, 
and the kinetic energy for the same reaction measured in the laboratory, we 


consider the deuterium-tritium fusion reaction 


2 4 1 
°D +:T —3He + hn 
2To be captured, the incident ion must approach the target nucleus to within a radius 


of 10~'°m, and must not have so much energy that it can escape the nuclear strong force 


existing at that radius. 
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in L, where D and T are the reacting deuterium and tritium ions, respectively; 
He and n are the helium nucleus and neutron, reaction products, respectively. 
The superscripts refer to the particle mass number, and are given in atomic 
mass units,? whereas the subscripts refer to the particle atomic number, or 
the number of protons in its nucleus. Here, the deuterium ion is accelerated 
in a positive direction towards the stationary (relative to the laboratory) tri- 


tium ion; upon interaction, the fusion reaction takes place. 


In the L frame, the total momentum of the the reacting particles is given 
by 
P=™pvp, 


with the total kinetic energy being 


1 
E= 5DUp: (A.1) 


where mp and vp are the deuterium mass and velocity measured in the L 
frame. In the CM frame, the total momentum of the reacting particles is 


correspondingly 


p=™MpUp+My0t7, 


=0, (A.2) 


with the total kinetic energy being 


a 
ES 5 (mpU_p + mrVy) , (A.3) 


where mr is the mass of the tritium ion, and where the over-line denotes 


quantities measured in the CM frame (U7 is the velocity of the tritium ion 


3An atomic mass unit is effectively the mass of a nucleon, or nuclear particle. Hence a 


particle with a given mass of 3 say, effectively consists of 3 nucleons. 
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in the CM frame, for example). Since the CM frame of reference moves in a 
positive direction relative to the laboratory, we denote the velocity of it, as 
measured in L, by vey. Furthermore, since only the deuterium ion is moving 


in L, then 


Up = Up — UGM (A.4) 


UT = —UcM, (A.5) 
where upon substitution of (A.4) and (A.5) into (A.2), we deduce that 


MD 
vem = | —— ] UD; 
MptmMmp 


consequently 


mr 
= | ——- A.6 
(<A) oe, (A.6) 
from (A.4). The total kinetic energy FE in the CM frame is then 
— MP 1 9 
Tae (een ee 
(— oars —) 5 PED 
= (—_) E (A.7) 
Mp + Mr 


from (A.1) and (A.3) - (A.6). 


A.2 Reaction Cross-Section in the CM and L 


Frames of Reference. 


Given (A.7), if we know the function o(£) representing the D-T reaction 


cross-section over a range of kinetic energies in the CM frame of reference, we 


204 


can scale this function’s dependent variable (£) by an appropriate parameter 
Az, where 
MD + mrp 


Gg] ——$——; 


MPL 


to give the equivalent function expressed in the laboratory (i.e. o(£) = 
a(A,E)). Figure A.1 shows the effect of scaling the CM generic D-T reaction 
cross-section curve (in black) to give both the D-T reaction (in blue) in the 
laboratory (where deuterium is incident upon tritium), and the T-D reaction 


(in red) in the laboratory (where tritium is incident upon deuterium). 


Reaction Cross-Sections for the D-T (T-D) Reaction 


Generic CM Cross Section Curve 
——— D-T Laboratory Frame Cross Section 
54 —— T-D Laboratory Frame Cross Section 


=(E) (barns) 
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Figure A.1: Reaction cross-section as a function of energy for the D-T reaction 
in the CM frame of reference (black), with the scaled cross-sections for both the 
D-T (red) and T-D (blue) in the laboratory frame of reference also shown. The 
cross-section peaks occur at E ~ 66 KeV for the CM frame, E =~ 110 KeV for 
D-T in the laboratory frame and E =~ 165 KeV for T-D in the laboratory frame. 


The original CM curve was taken from [4] and [5]. 
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Appendix B 


Magnetic Field Effects 


In a neutron tube no external magnetic field is applied during its operation, 
and so the only possible contribution to B, on the right hand side of (2.13a), 
is the magnetic field induced by the current of charged particles themselves 
as they flow across the region. Such a magnetic field is naturally induced by 
the electric field advection, caused by this passage of charged particles, and 
is a consequence of the continuity equation (2.13b) and of the electric field 


divergence, (2.13c); it can be derived as follows. 


e Magnetic field induction due to current flow 


By equating the time derivative of (2.13c) to (2.13b), we have 


OE 
Vv : {ort eos} = 0, 


so that upon integration we obtain, 


OE 
pv + co" = alt), (B.1) 


where the vector function g(t) can be a function of time only. By 


comparing (B.1) to the well known differential form of Ampére’s law 
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(noting the definition (2.16)), 
OE 
Lo (3 + a) =VxB (B.2) 


(see [14], for example), we see that g(t) = aa WV XB, with the implication 
that V x B must also be a function of time only; here, jo is a scaling 
parameter known as the permeability of free space. Re-writing (B.2) 


using (2.13c) we obtain 
E 
(vv -E+ =) =¢V x B; (B:3) 


where c = /1/p0E0 is the speed of light in vacuum, with the immediate 
conclusion being that the flow of charged particles naturally introduces 
electric field advection, which in turn introduces a magnetic field that 
is perpendicular to the particle flow. The curl of this field is purely 
time dependent and by writing f(t) = g(t)/e9 = c?V x B, we can say 


in two Cartesian dimensions that 


f(t) = filt)e+ folb)y, 
=?VxB 


OB OB OB OB 
a Be 2 Vy \s _ Zl z\ a 
=" ( Oy Oz ji ( Ox Oz )4 ; ee) 


In this thesis, we are only concerned with planar tube models in one 


and two dimensions, with a planar, two-dimensional (in the indepen- 
dent variables x and y) region being effectively the cross-section of a 
prism extending infinitely into the z direction. Hence, in this frame- 


work, partial derivatives in the z direction are zero, implying that the 
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functions f;(t) and fo(t) are given by 


—— B. 
lt) = 0, (Ba) 
OB 
Be Bee B.5b 
f(t) = 2, (Bb) 
since 0B, /0z = 0B,/0z = 0. We can then say that 
OB, =yfi(t) + Fi, (B.6) 
from (B.5a), and 
cB, a —2 fo(t) + Fo, (B.7) 


from (B.5b), where -¥; and .¥9 are arbitrary functions of (x,t) and 
(y,t), respectively. By combining (B.6) and (B.7), the z component of 


the self-induced magnetic field can be written 


B.=4(vhlt) —2hl) +90), (B8) 
where Y is an arbitrary function of time, clearly showing that the elec- 
tric field advection induces the perpendicular magnetic field component 
B, = B,(x,y). This component is small, indicated by the presence of 
1/c? on the right hand side of (B.8). Moreover, the presence of c? on 
the right hand side of (B.3) offsets the small magnitude of the self- 
induced field B, making the vector function f(t) comparable to the 
electric field advection, thus balancing (B.3). The implication here is 
that both V x B and B, are small functions (the presence of c? in 
both cases clearly indicates this) relative to the electric field advection 


propagated by the charged particle flow. 


We now seek to futher show that the effects of the small self-induced 
magnetic field on charged particle flow are also negligible in comparison 


to the effects caused by the application of the external electric field. 
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e Effect of self-induced magnetic field on particle acceleration 


The component of accelerative force acting on the charged particle 


beam as a consequence of any magnetic fields is given by 
F,=qv xB, (B.9) 


from (2.13a), or in this case, because the only magnetic field present is 


that due to self-induction, 


since the self-induced magnetic field only has a component perpendic- 
ular to the direction of particle motion (given by (B.8)), and where v, 


and v, are the x and y components of the particle velocity, respectively. 
One-dimensional case 


In one dimension, acceleration in the 7 direction can be ignored, and 
additionally, there is no component of particle velocity in this direction. 
This implies (from (B.10)) that in one dimension, acceleration due to 


self-induced magnetic field can be ignored. 
Two-dimensional case 


From (B.10), it can be seen that in two dimensions, the magnitude 
of the self-induced magnetic field can cause an accelerative force in 
the x—y plane. However, as a consequence of the presence of 1/c? in 
(B.8), we know that the self-induced magnetic field component, B., 
must be small in comparison to any electric field advection caused by 
the charged particle flow, but we do not know the sizes of the functions 


fi(t), fo(t) and Y(t). Nonetheless, an estimate of the magnitude of the 
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self-induced magnetic field can be determined from the Biot-Savart law 
({14], p296) where its magnitude some distance r radially away from a 
long, current-carrying conductor (such as a neutron tube acceleration 
gap) carrying a current J, is given by 

|B] = ae (B.11) 
From (B.9), the magnitude of the component of the force acting on the 
beam of charged particles, due to the self-induced magnetic field at r 


is then 


[Fs] = |B] |v| 


f 
— qbol |v) (B.12) 
2ar 


Choosing r to be of the order of the ion beam radius (r ~ 107 m), we 
find that for a typical 1A ion current, with the mean ion velocity being 
that of a particle with 50KeV energy (the typical mean energy of a 
deuteron within a neutron tube acceleration region), the magnitude of 
the force due to the self-induced magnetic field at one ion beam radius 
is |Fp| ~ 7 x 107!" N. In contrast, the magnitude of the component 
of the force acting on the charged particles due to the applied electric 
field, 
[Fx| = q|B| 


(from (2.13a)), is some five orders of magnitude greater for a typical 
acceleration potential difference of ~ 100 kV and acceleration gap of 
~10-?m. With this in mind, the contribution to particle acceleration 


in (2.13a) (perpendicular to the direction of particle flow) due to the 
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self-induced magnetic field is deemed negligible and, in conjunction to 


experimental experience, allows it to be disregarded in this study. 


Time varying magnetic field influence. 


The link between the electric and magnetic fields via the time derivative 


of the magnetic vector potential in (2.13e) implies that 


dB 
VxE=-3, (B.13) 


from (2.13f). Consequently, in the presence of a time varying magnetic 
field (as is the case with charged particle flow, shown in (B.8)), the 
electric field defined by (2.13f) clearly has the non-zero curl, (B.13), and 
cannot truly be considered a conservative field. However, by examining 
(B.13), we can show that no matter how large the rate of change of 
magnetic field at the wavefront is, it can be neglected, validating a 
conservative field approximation. We do this by expanding the curl 
term on the left hand side of (B.13), giving 


_ (OE. OE,\  , (OE, OE:\ . + (OE, OE; 
vxE=i(5 =) +3(F ma) +k (SS =) 


(B.14) 


where FE, represents the component of E in the x direction etc. Given 
that E(x, t) is, at most, a two-dimensional field (in this thesis), its z 
component does not exist; moreover, derivatives in the z direction are 


zero. With these things in mind, (B.14) can be reduced to 
~ { OF, Ey 
VxE=k(52-5 ). (B.15) 


showing that V x E, and consequently 0OB/Ot (consistently with (B.8)) 


only have a non-zero component pointing out of the x—y plane. There- 
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fore, in the x—y plane we must conclude 
Vx E=}0, 


from (B.15), implying that the electric field is conservative (with no 
vorticity) and that 
E=-—V¢ (B.16) 


must be true for some ¢, from (2.13e). Hence, for the purposes of this 
study, (2.18) is considered the definition the electric field E(x, t) for 


the scalar potential $(x, t). 
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